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Preface to the Second Edition 



Since the first edition of this book was published, there have been major improve- 
ments in symbolic mathematical languages such as Maple™ and Mathematica™ 
and this has opened up the possibility of solving considerably more complex and hence 
interesting and realistic elasticity problems as classroom examples. It also enables the 
student to focus on the formulation of the problem (e.g. the appropriate governing 
equations and boundary conditions) rather than on the algebraic manipulations, with 
a consequent improvement in insight into the subject and in motivation. During the 
past 10 years I have developed files in Maple and Mathematica to facilitate this pro- 
cess, notably electronic versions of the Tables in the present Chapters 19 and 20 and 
of the recurrence relations for generating spherical harmonics. One purpose of this 
new edition is to make this electronic material available to the reader through the 
Kluwer website www.elasticity.org. I hope that readers will make use of this resource 
and report back to me any aspects of the electronic material that could benefit from 
improvement or extension. Some hints about the use of this material are contained in 
Appendix A. Those who have never used Maple or Mathematica will find that it takes 
only a few hours of trial and error to learn how to write programs to solve boundary 
value problems in elasticity. 

I have also taken the opportunity to include substantially more material in the 
second edition — notably three chapters on antiplane stress systems, including Saint- 
Venant torsion and bending and an expanded section on three-dimensional problems 
in spherical and cylindrical coordinate systems, including axisymmetric torsion of 
bars of non-uniform circular cross-section. 

Finally, I have greatly expanded the number of end-of-chapter problems. Some of 
these problems are quite challenging, indeed several were the subject of substantial 
technical papers within the not too distant past, but they can all be solved in a few 
hours using Maple or Mathematica. A full set of solutions to these problems is in 
preparation and will be made available to bona fide instructors on request. 

J. R. Barber 
Ann Arbor 
2002 




Preface to the First Edition 



The subject of Elasticity can be approached from several points of view, depending 
on whether the practitioner is principally interested in the mathematical structure of 
the subject or in its use in engineering applications and in the latter case, whether 
essentially numerical or analytical methods are envisaged as the solution method. My 
first introduction to the subject was in response to a need for information about a 
specific problem in Tribology. As a practising Engineer with a background only in 
elementary Strength of Materials, I approached that problem initially using the con- 
cepts of concentrated forces and superposition. Today, with a rather more extensive 
knowledge of analytical techniques in Elasticity, I still find it helpful to go back to 
these roots in the elementary theory and think through a problem physically as well 
as mathematically, whenever some new and unexpected feature presents difficulties in 
research. This way of thinking will be found to permeate this book. My engineering 
background will also reveal itself in a tendency to work examples through to final 
expressions for stresses and displacements, rather than leave the derivation at a point 
where the remaining manipulations would be routine. 

With the practical engineering reader in mind, I have endeavoured to keep to 
a minimum any dependence on previous knowledge of Solid Mechanics, Continuum 
Mechanics or Mathematics. Most of the text should be readily intelligible to a reader 
with an undergraduate background of one or two courses in elementary Strength of 
Materials and a rudimentary knowledge of partial differentiation. Cartesian tensor 
notation and the summation convention are used in a few places to shorten the deriva- 
tion of some general results, but these sections are carefully explained, so as to be 
self-explanatory. 

The book is based on a one semester graduate course on Linear Elasticity that I 
have taught at the University of Michigan since 1983. In such a restricted format, it 
is clearly necessary to make some difficult choices about which topics to include and, 
more significantly, which to exclude. The most significant exclusion is the classical 
complex variable solution of two-dimensional Elasticity, which, if it were to be ade- 
quately treated including essential mathematical preliminaries, would need most of a 
book of this length to itself. Instead, I have chosen to restrict the two-dimensional 
treatment to the more traditional real stress function approach, so as to leave room for 
a substantial amount of material on three-dimensional problems, which are arguably 
closer to the frontier of current research. 
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Modern practitioners of Elasticity are necessarily influenced by developments in 
numerical methods, which promise to solve all problems with no more information 
about the subject than is needed to formulate the description of a representative 
element of material in a relatively simple state of stress. As a researcher in Solid 
Mechanics, with a primary interest in the physical behaviour of the systems I am 
investigating, rather than in the mathematical structure of the solutions, I have fre- 
quently had recourse to numerical methods of all types and have tended to adopt 
the pragmatic criterion that the best method is that which gives the most convincing 
and accurate result in the shortest time. In this context, ‘convincing’ means that the 
solution should be capable of being checked against reliable closed-form solutions in 
suitable limiting cases and that it is demonstrably stable and in some sense conver- 
gent. Measured against these criteria, the ‘best’ solution to many practical problems 
is often not a direct numerical method, such as the finite element method, but rather 
one involving some significant analytical steps before the final numerical evaluation. 
This is particularly true in three-dimensional problems, where direct numerical meth- 
ods are extremely computer-intensive if any reasonably accuracy is required, and in 
problems involving infinite or semi-infinite domains, discontinuities, bonded or con- 
tacting material interfaces or theoretically singular stress fields. By contrast, I would 
immediately opt for a finite element solution of any two-dimensional problem involv- 
ing finite bodies with relatively smooth contours, unless it happened to fall into the 
(surprisingly wide) class of problems to which the solution can be written down in 
closed form. The reader will therefore find my choice of topics significantly biassed 
towards those fields identified above where analytical methods are most useful. 

I have provided a representative selection of problems suitable for class use at the 
end of most of the chapters. Many texts on Elasticity contain problems which offer a 
candidate stress function and invite the student to ‘verify’ that it defines the solution 
to a given problem. Students invariably raise the question ‘How would we know to 
choose that form if we were not given it in advance?’ I have tried wherever possible to 
avoid this by expressing the problems in the form they would arise in Engineering — 
i.e. as a body of a given geometry subjected to prescribed loading. This in turn has 
required me to write the text in such a way that the student can approach problems 
deductively. I have also generally opted for explaining difficulties that might arise in 
an ‘obvious’ approach to the problem, rather than steering the reader around them 
in the interests of brevity. 

Even conceptually straightforward problems in three-dimensional Elasticity tend 
to be algebraically complicated, because of the dependence of the results on Poisson’s 
ratio, which is absent in two-dimensions, except where body forces etc. are involved. 
With the greater general availablity of symbolic processors, this will soon cease to 
be a difficulty, but those without access to such methods might like to simplify some 
of the three-dimensional problems by restricting the solution to the case v = 0. At 
the University of Michigan, we are fortunate in having a very sophisticated computer 
environment, so I can ask my students to carry the problems as far as plotting graphs 
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of stress distributions and performing parametric studies to determine optimum de- 
signs for components. From the engineering perspective, this leads to a significantly 
greater insight into the nature of the subject, as well as serving as an educational 
motivation for those students outside the classical Mechanics field. 

Over the years, the students attending my course have been drawn from depart- 
ments of Mechanical, Civil, Aerospace and Materials Engineering as well as Naval 
Architecture, Applied Mechanics and Mathematics and have therefore generally rep- 
resented quite a diverse set of points of view. The questions raised by this very 
talented group of people have often resulted in vigorous classroom discussions which 
have had a major impact on my perspective on the subject. It is therefore to my 
graduate students, past and present, that this book is dedicated. 

It is traditional at this point also to thank a multitude of colleagues, assistants 
and family members for their contributions and support. However, since I developed 
the text as part of my normal teaching duties and typset it myself using IATj^X, many 
of these credits would be out of place. My wife, Maria Comninou, may disagree, 
but as far as I can judge, my temperament was not significantly soured during the 
gestation period, but she does deserve my thanks for her technical input during times 
that we have collaborated on research problems. 

I am particularly grateful to the Rector and Fellows of Lincoln College, Oxford, 
who provided me with a stimulating scholarly environment during an important stage 
in the writing, and to David Hills of that College who made my visit there possible. 

The line drawings were prepared by Rodney Hill and Figures 13.2, 13.3 are repro- 
duced by kind permission of the Royal Society of London. I should also like to thank 
John Dundurs for permission to use Table 9.1, whose usefulness has been proved by 
his students at Northwestern University and many others in the Elasticity community. 

J.R. Barber 
Ann Arbor 
1991 
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GENERAL CONSIDERATIONS 




Chapter 1 

INTRODUCTION 



The subject of Elasticity is concerned with the determination of the stresses and 
displacements in a body as a result of applied mechanical or thermal loads, for those 
cases in which the body reverts to its original state on the removal of the loads. 
In this book, we shall further restrict attention to the case of linear infinitesimal 
elasticity, in which the stresses and displacements are linearly proportional to the 
applied loads and the displacements are small in comparison with the characteristic 
length dimensions of the body. These restrictions ensure that linear superposition 
can be used and enable us to employ a wide range of series and transform techniques 
which are not available for non-linear problems. 

Most engineers first encounter problems of this kind in the context of the subject 
known as Mechanics of Materials, which is an important constituent of most under- 
graduate engineering curricula. Mechanics of Materials differs from Elasticity in that 
various plausible but unsubstantiated assumptions are made about the deformation 
process in the course of the analysis. A typical example is the assumption that plane 
sections remain plane in the bending of a slender beam. Elasticity makes no such 
assumptions, but attempts to develop the solution directly and rigorously from its 
first principles, which are Newton’s laws of motion, Euclidian geometry and Hooke’s 
law. Approximations are often introduced towards the end of the solution, but these 
are mathematical approximations used to obtain solutions of the governing equations 
rather than physical approximations that impose artificial and strictly unjustifiable 
constraints on the permissible deformation field. 

However, it would be a mistake to draw too firm a distinction between the two 
approaches, since practitioners of each have much to learn from the other. Mechanics 
of Materials, with its emphasis on physical reasoning and a full exploration of the 
practical consequences of the results, is often able to provide insights into the problem 
that are less easily obtained from a purely mathematical perspective. Indeed, we 
shall make extensive use of physical parallels in this book and pursue many problems 
to conclusions relevant to practical applications, with the hope of deepening the 
reader’s understanding of the underlying structure of the subject. Conversely, the 
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CHAPTER 1. INTRODUCTION 



mathematical rigour of Elasticity gives us greater confidence in the results, since, 
even when we have to resort to an approximate solution, we can usually estimate 
its accuracy with some confidence — something that is very difficult to do with the 
physical approximations used in Mechanics of Materials 1 . Also, there is little to be 
said for using an ad hoc approach when, as is often the case, a more rigorous treatment 
presents no serious difficulty. 



1.1 Notation for stress and displacement 

It is assumed that the reader is more or less familiar with the concept of stress and 
strain from elementary courses on Mechanics of Materials. This section is intended 
to introduce the notation used, to refresh the reader’s memory about some important 
ideas, and to record some elementary but useful results. 

1.1.1 Stress 

Components of stress will all be denoted by the symbol o with appropriate suffices. 
The second suffix denotes the direction of the stress component and the first the 
direction of the outward normal to the surface upon which it acts. This notation is 
illustrated in Figure 1.1 for the Cartesian coordinate system x,y,z. 




Figure 1.1: Notation for stress components. 

Notice that one consequence of this notation is that normal (i.e. tensile and 
compressive) stresses have both suffices the same (e.g. o xx ,o yy ,o Z2 in Figure 1.1) and 
are positive when tensile. The remaining six stress components in Figure 1.1 (i.e. 
o xy , CTyi, Qzy, ^zxi &xz) have two different suffices and are shear stresses. 

In fact, the only practical way to examine the effect of these approximations is to relax them, 
by considering the same problem, or maybe a simpler problem with similar features, in the context 
of the theory of Elasticity. 
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Books on Mechanics of Materials often use the symbol r for shear stress, whilst 
retaining a for normal stress. However, there is no need for a different symbol, since 
the suffices enable us to distinguish normal from shear stress components. Also, we 
shall find that the use of a single symbol with appropriate suffices permits matrix 
methods to be used in many derivations and introduces considerable economies in 
the notation for general results. 

The equilibrium of moments acting on the block in Figure 1.1 requires that 

(Jxy — &yx i rTy ^ — U zy &nd Ozx ™ @xz • (1.1) 

This has the incidental advantage of rendering mistakes about the order of suffices 
harmless! (In fact, a few books use the opposite convention.) Readers who have not 
encountered three-dimensional problems before should note that there are two shear 
stress components on each surface and one normal stress component. There are some 
circumstances in which it is convenient to combine the two shear stresses on a given 
plane into a two-dimensional vector in the plane — i.e. to refer to the resultant shear 
stress on the plane. An elementary situation where this is helpful is in the Mechanics 
of Materials problem of determining the distribution of shear stress on the cross- 
section of a beam due to a transverse shear force 2 . For example, we note that in this 
case, the resultant shear stress on the plane must be tangential to the edge at the 
boundary of the cross section, since the shear stress complementary to the component 
normal to the edge acts on the traction-free surface of the beam and must therefore 
be zero. This of course is why the shear stress in a thin-walled section tends to follow 
the direction of the wall. 

We shall refer to a plane normal to the x-direction as an ‘x-plane’ etc. The only 
stress components which act on an x-plane are those which have an x as the first 
suffix (This is an immediate consequence of the definition). 

Notice also that any x-plane can be defined by an equation of the form x = C, 
where C is a constant. More precisely, we can define a ‘positive x-plane’ as a plane for 
which the positive x-direction is the outward normal and a ‘negative x-plane’ as one 
for which it is the inward normal. This distinction can be expressed mathematically 
in terms of inequalities. Thus, if part of the boundary of a solid is expressible as 
X = C, the solid must locally occupy one side or the other of this plane. If the domain 
of the solid is locally described by x < C, the bounding surface is a positive x-plane, 
whereas if it is described by x > C } the bounding surface is a negative x-plane. 

The discussion in this section suggests a useful formalism for correctly defining the 
boundary conditions in problems where the boundaries are parallel to the coordinate 
axes. We first identify the equations which define the boundaries of the solid and then 
write down the three traction components which act on each boundary. For example, 
suppose we have a rectangular solid defined by the inequalities 0<x<a, 0 < y < 
b, 0 < z < c. It is clear that the surface y = 6 is a positive y-plane and we deduce 

2 For the Elasticity solution of this problem, see Chapter 17. 
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immediately that the corresponding traction boundary conditions will involve the 
stress components a yx ,Oyy,o yz — i.e. the three components that have y as the first 
suffix. This procedure insures against the common student mistake of assuming (for 
example) that the component a xx must be zero if the surface y = b is to be traction 
free. (Note : Don’t assume that this mistake is too obvious for you to fall into it. 
When the problem is geometrically or algebraically very complicated, it is only too 
easy to get distracted.) 

Stress components can be defined in the same way for other systems of orthogonal 
coordinates. For example, components for the system of cylindrical polar coordinates 
(r,9,z) are shown in Figure 1.2. 




Figure 1.2: Stress components in polar coordinates. 

(This is a case where the definition of the ‘0-plane’ through an equation, 9 ~ C, is 
easier to comprehend than ‘the plane normal to the 0-direction’. However, note that 
the 0-direction is the direction in which a particle would move if 0 were increased 
with r, z constant.) 

1.1.2 Index and vector notation 

Many authors use the notation x t , i 2 ,l 3 , in place of I, y, z for the Cartesian coordinate 
system, in which case the stress components are written 011 , 012 , etc. (<Tnxi> o XlX2 
etc. would obviously be too cumbersome). This notation has the particular advantage 
that in combination with the ‘summation convention’ it permits general results to be 
written and manipulated in a concise and elegant form. The summation convention 
lays down that any term in which the same latin suffix occurs twice stands for the 
sum of all the terms obtained by giving this suffix each of its possible values. For 
example <7,i is interpreted as 

3 

0ji =■ = On + 022 + 033 

t=l 



( 1 . 2 ) 
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and 



dui _ du\ du 2 du 3 

dxi dxi dx 2 dx 3 



(1.3) 



On the other hand, some results are more conveniently manipulated or expressed 
in vector notation, for which we define the position vector 



R = ix + jy + kz , 

where i, j, k are unit vectors in directions x,y,z respectively. 
Important results will be given in each form wherever possible. 



(1.4) 



1.1.3 Vectors, tensors and transformation rules 

Vectors can be conceived in a mathematical sense as ordered sets of numbers or 
in a physical sense as mathematical representations of quantities characterized by 
magnitude and direction. A link between these concepts is provided by the transfor- 
mation rules. Suppose we know the components ( u x ,u y ) of the vector it in a given 
two-dimensional Cartesian coordinate system (x, y) and we wish to determine the 
components ( u ' x , u' y ) in a new system ( x y') which is inclined to ( x , y) at an angle 6 
in the counterclockwise direction as shown in Figure 1.3. The required components 
are 



u' x — u x cos 0 T Uy sin 6 (1.5) 

u' y = u y cos 6 — u x sin 0 . (1.6) 




Figure 1.3: The coordinate systems x, y and x',y'. 

We could define a vector as an entity, described by its components in a specified 
Cartesian coordinate system, which transforms into other coordinate systems accord- 
ing to rules like equations (1.5, 1.6) — i.e. as an ordered set of numbers which obey 
the transformation rules (1.5, 1.6). The idea of magnitude and direction could then 
be introduced by noting that we can always choose 0 such that (i) u' y = 0 and (ii) 
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u' x > 0. The corresponding direction x' is then the direction of the resultant vector 
and the component u' x is its magnitude. 

Now stresses have two suffices and components associated with all possible com- 
binations of two coordinate directions, though we note that equation (1.1) shows that 
the order of the suffices is immaterial. (Another way of stating this is that the matrix 
of stress components a ; ; is always symmetric). The stress components satisfy a more 
complicated set of transformation rules which in the two-dimensional case are those 
associated with Mohr's circle — i.e. 



Gx'x' 


— o xx cos 2 6 + o yy sin 2 0 + 2 a xy sin 0 cos 6 


(1.7) 


Ox‘y' 


= o xy (c os 2 0 - sin 2 6) + (o yy — o xx ) sin 9 cos 6 


(1.8) 


a y'y' 


= o yy cos 2 9 + o xx sin 2 9 — 2 o xy sin 6 cos 6 . 


(1.9) 



As in the case of vectors we can define a mathematical entity which has a matrix 
of components in any given Cartesian coordinate system and which transforms into 
other such coordinate systems according to rules like (1.7-1 .9). Such quantities are 
called second order Cartesian tensors. 

We know from Mohr's circle that we can always choose 6 such that o x 'y> =0, in 
which case the directions x\ y' are referred to as principal directions and the compo- 
nents ay , Oyty' as principal stresses. Thus another way to characterize a second order 
Cartesian tensor is as a quantity defined by a set of orthogonal principal directions 
and a corresponding set of principal values. 

As with vectors, a pragmatic motivation for abstracting the mathematical prop- 
erties from the physical quantities which exhibit them is that many different physical 
quantities are naturally represented as second order Cartesian tensors. Apart from 
stress and strain, some commonly occurring examples are the second moments of area 
of a beam cross section ( I xx , 7^, I yy ), the second partial derivatives of a scalar func- 
tion (d 2 f/dx 2 \ d 2 fjdxdy\ d 2 f /dy 2 ) and the influence coefficient matrix Cij defining 
the displacement u due to a force F for a linear elastic system, i.e. 

u-i — CijFj , (1.10) 

where the summation convention is implied. 

It is a fairly straightforward matter to prove that each of these quantities obeys 
transformation rules like (1.7-1. 9). It follows immediately (for example) that every 
beam cross section has two orthogonal principal axes of bending about which the two 
principal second moments are respectively the maximum and minimum for the cross 
section. 

A special tensor of some interest is that for which the Mohr's circle degenerates 
to a point. In the case of stresses, this corresponds to a state of hydrostatic stress, 
so-called because a fluid at rest cannot support shear stress (the constitutive law for 
a fluid relates velocity gradient to shear stress) and hence o x < y < — 0 for all 0. The 
only Mohr's circle which satisfies this condition is one of zero radius, from which we 
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deduce immediately that all directions are principal directions and that the principal 
values are all equal. In the case of the fluid, we obtain the well-known result that the 
pressure in a fluid at rest is equal in all directions. 

It is instructive to consider this result in the context of other systems involving 
tensors. For example, consider the second moments of area for the square cross section 
shown in Figure 1.4. By symmetry, we know that Ox, Oy are principal directions and 
that the two principal second moments are both equal to a‘*/12. It follows immediately 
that the Mohr’s circle has zero radius and hence that the second moment about any 
other axis must also be a 4 / 12 — a result which is not obvious from an examination 
of the section. 



t r 



a 










a 


O 







Figure 1.4: A beam of square cross-section. 

As a second example, Figure 1.5 shows an elastic system consisting of three iden- 
tical but arbitrary structures connecting a point, P, to a rigid support, the structures 
being inclined to each other at angles of 120°. 

The structures each have elastic properties expressible in the form of an influence 
function matrix as in equation (1.10) and are generally such that the displacement 
u is not colinear with the force F. However, the overall influence function matrix 
for the system has the same properties in three different coordinate systems inclined 
to each other at 120°, since a rotation of the Figure through 120° leaves the system 
unchanged. The only Mohr’s circle which gives equal components after a rotation of 
120° is that of zero radius. We therefore conclude that the support system of Figure 
1.5 is such that (i) the displacement of P always has the same direction as the force 
F and (ii) the stiffness or compliance of the system is the same in all directions 3 . 

3 A similar argument can be used to show that if a laminated fibre-reinforced composite is laid 
up with equal numbers of identical, but not necessarily symmetrical, laminas in each of 3 or more 
equispaced orientations, it must be elastically isotropic within the plane. This proof depends on the 
properties of th e fourth order Cartesian tensor Cy*; describing the stress-strain relation e,j ~CijklOkl 
of the laminas. 
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Figure 1.5: Support structure with 3 similar but unsymmetrical components. 

These two examples illustrate that there is sometimes an advantage to be gained 
from considering a disparate physical problem that shares a common mathematical 
structure with that under investigation. 

1.1.4 Principal stresses and Von Mises stress 

One of the principal reasons for performing elasticity calculations is to determine 
when and where an engineering component will fail. Theories of material failure 
are beyond the scope of this book, but the most widely used criteria are the Von 
Mises distortion energy criterion for ductile materials and the maximum tensile stress 
criterion for brittle materials 4 . 

Brittle materials typically fracture when the maximum tensile stress reaches a 
critical value and hence we need to be able to calculate the maximum principal stress 
o i . For two dimensional problems, the principal stresses can be found by using the 
condition a x > v < — 0 in equation (1.8) to determine the inclination 8 of the principal 
directions. Substituting into (1.7, 1.9) then yields the well known results 

w = ± + < ■ (mi 

In most problems, the maximum stresses occur at the boundaries where shear 
tractions are usually zero. Thus, even in three-dimensional problems, the determi- 
nation of the maximum tensile stress often involves only a two-dimensional stress 

4 J.R. Barber, Intermediate Mechanics of Materials, McGraw-Hill, New York (2000), §2.2. 
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transformation. However, the principal stresses are easily obtained in the fully three- 
dimensional case from the results 



CTi = 


± + 2 ( 
3 3V 


CO 

1 

<M_ 


j COS (j) 




h 2 / 


EE, TT ' 


\ f , 271 


cr 2 = 


3 + 3 ( 


\jl\ - 3 I 2 


) COS ((» + — 




h 2 / 


EE TV ' 




a 3 = 


"3 + 3 ( 


\J If - 3 I 2 


) cos (d+Y 



( 1 . 12 ) 

(1.13) 

(1.14) 



where 



and 



• 1 ( 2/f - 9IJ 2 + 27/ 3 

0 = - arccos ^- ( - > F - 3 - 3?r 



(1.15) 



h 

h 

h 



= a. 



xx Qyy ^ ^ z 



®xx®yy "t" GyxjGzz "t" Ozz^xx & ™ ^nr ® 



= < 7 . 



'xy 

__ 2 _ 2 __ 
xx^yyOzz &xx&yz ^yy^zx ®zz® 



yz v zx 
xy 20xy CT yz ( 7zx 



(1.16) 

(1.17) 

(1.18) 



The quantities 7), I 2 , 13 are known as stress invariants because for a given stress state 
they are the same in all coordinate systems. 

Von Mises theory states that a ductile material will yield when the strain energy 
of distortion per unit volume reaches a certain critical value. This leads to the yield 
criterion 

o E = \Jly~ 3/2 

= \J olx + crly + <J 2 ZZ - OxxOyy ~ Oyy° Z z ~ Ozz°xx + 3 0 2 xy + 3 + 3(jj x = Sy , ( 1 . 1 9) 

where Sy is the yield stress in uniaxial tension. The quantity a E is known as the 
equivalent tensile stress or the Von Mises stress. 

The Maple and Mathematica files ‘principalstresses’ use equations (1.12-1.19) to 
calculate the principal stresses and the Von Mises stress from a given set of stress 
components. 



1.1.5 Displacement 

The displacement of a particle P is a vector u representing the difference between 
the final and the initial position of P — i.e. it is the distance which P moves during 
the deformation. The components of u are denoted by appropriate suffices — e.g. 
u x , u v , u z , so that 

u = iu x + ju y + ku z . (1.20) 

The deformation of a body is completely defined if we know the displacement of 
its every particle. Notice however that there is a class of displacements which do 
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not involve deformation — the so-called ‘rigid-body displacements’. A typical case is 
where all the particles of the body have the same displacement. The name arises, of 
course, because rigid-body displacement is the only class of displacement that can be 
experienced by a rigid body. 



1.2 Strains and their relation to displacements 

Components of strain will be denoted by the symbol, e, with appropriate suffices 
(e.g. e xx ,e xy ). As in the case of stress, no special symbol is required for shear strain, 
though we shall see below that the quantity defined in most elementary texts (and 
usually denoted by 7 ) differs from that used in the mathematical theory of Elasticity 
by a factor of 2. A major advantage of this definition is that it makes the strain, 
e, a second order Cartesian Tensor (see §1.1.3 above). We shall demonstrate this by 
establishing transformation rules for strain similar to equations (1.7- 1.9) in §1.2.4 
below. 

1.2.1 Tensile strain 

Students usually first encounter the concept of strain in elementary Mechanics of 
Materials as the ratio of extension to original length and are sometimes confused by 
the apparently totally different definition used in more mathematical treatments of 
solid mechanics. We shall discuss here the connection between the two definitions — 
partly for completeness, and partly because the physical insight that can be developed 
in the simple problems of Mechanics of Materials is very useful if it can be carried 
over into more difficult problems. 

Figure 1.6 shows a bar of original length L and density p hanging from the ceiling. 
Suppose we are asked to find how much it increases in length under the loading of its 
own weight. 




Figure 1.6: The bar suspended from the ceiling. 
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It is easily shown that the tensile stress a xx at the point P, distance x from the 
ceiling, is 



o xx = pg(L - x) , 

where g is the acceleration due to gravity, and hence from Hooke’s law, 

pg{L - x ) 

e xx — c 



( 1 . 21 ) 



( 1 . 22 ) 



However, the strain varies continuously over the length of the bar and hence we 
can only apply the Mechanics of Materials definition if we examine an infinitesimal 
piece of the bar over which the strain can be regarded as sensibly constant. 

We describe the deformation in terms of the downward displacement u x which 
depends upon x and consider that part of the bar between x and x+6x, denoted by 
PQ in Figure 1.7. 



x 



x+8x 



P 




Figure 1.7: Infinitesimal section of the bar. 



After the deformation, PQ must have extended by u x (Q) — u x (P) and hence the local 
value of ‘Mechanics of Materials' tensile strain is 



u x (Q) - u x (P) _ u x (x + 5x) - u x (x ) 
Sx Sx 



(1.23) 



Taking the limit as 5x — > 0» we obtain the definition 

_ du x 

&XX Q 

OX 



(1.24) 



Corresponding definitions can be developed in three-dimensional problems for the 
other normal strain components, i.e. 

_ du y _ du z 

eyy ~~d^ ' e “ “ !h • 



(1.25) 
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Notice how easy the problem of Figure 1.6 becomes when we use these definitions. 
We get 



du x _ pg[L - x) 
dx E 



(1.26) 



from (1.22, 1.24) and hence 



u x = 



pq(2Lx - x 2 ) 

"2E 1+A 



(1.27) 



where A is an arbitrary constant of integration which expresses the fact that our 
knowledge of the stresses and hence the strains in the body is not sufficient to deter- 
mine its position in space. In fact, A represents an arbitrary rigid-body displacement. 
In this case we need to use the fact that the top of the bar is joined to a supposedly 
rigid ceiling — i.e. Uz(0) = 0 and hence A = 0 from (1.27). 



1.2.2 Rotation and shear strain 

Noting that the two z’s in e xx correspond to those in its definition du x /dx , it is natural 
to seek a connection between the shear strain e xy and one or both of the derivatives 
du x /dy,du y /dx. As a first step, we shall discuss the geometrical interpretation of 
these derivatives. 




Figure 1.8: Rotation of a line segment. 



Figure 1.8 shows a line segment PQ of length 6x, aligned with the a:-axis, the two 
ends of which are displaced in the y-direction. Clearly if u y (Q) / u y (P), the lineP<2 
will be rotated by these displacements and if the angle of rotation is small it can be 
written 



<t> = 



u y (x 4- Sx) - u y (x) 
6x 



( 1 . 28 ) 



(anticlockwise positive). 
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Proceeding to the limit as 5x — > 0, we have 



(1.29) 



Thus, du y /dx is the angle through which a line originally in the x-direction rotates 
towards the y-direction during the deformation 3 . 

Now, if PQ is a line drawn on the surface of an elastic solid, the occurence of 
a rotation (p does not necessarily indicate that the solid is deformed — we could 
rotate the line simply by rotating the solid as a rigid body. To investigate this matter 
further, we imagine drawing a series of lines at different angles through the point P 
as shown in Figure 1.9(a). 




(a) 





Figure 1.9: Rotation of lines at point P; (a) Original state, ( b ) Rigid-body 
rotation; (c) Rotation and deformation. 

If the vicinity of the point P suffers merely a local rigid-body rotation, all the 
lines will rotate through the same angle and retain the same relative inclinations 
as shown in 1.9 (b). However, if different lines rotate through different angles, as in 
1.9(c), the body must have been deformed. We shall show in the next section that 
the rotations of the lines in Figure 1.9(c) are not independent and a consideration of 
their interdependency leads naturally to a definition of shear strain. 

1.2.3 Transformation of coordinates 

Suppose we knew the displacement components u x ,u y throughout the body and 
wished to find the rotation, (p, of the line PQ in Figure 1.10, which is inclined at 
an angle 0 to the x-axis. 

5 Students of Mechanics of Materials will have already used a similar result when they express 
the slope of a beam as du/dx, where x is the distance along the axis of the beam and u is the 
transversedisplacement. 
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We construct a new axis system Ox'y' with Ox' parallel to PQ as shown, i 
case we can argue as above that PQ rotates anticlockwise through the angle 



HPQ) = 




Figure 1.10: Rotation of a line inclined at angle 9. 
Furthermore, we have 

■k ” V i ' 

“ ( i |x +3 '|) i ' = U S + ^'| i 

= cos0— - +siri0— , 
ox ay 



and by similar arguments, 



u x' = cos 9 4- u v sin 9 



u y > = u y cos 9 — u x sin 9 . 

Substituting these results into equation (1.30), we find 

<t>(PQ) = ^ cos 9 + sin 9-^ j ( u y cos 9 - u x sin 0) 



Oily 2 ri du x . 2 n 

cos 2 9 - — sin 2 0 4- 

ox 



dy 



du y _ dux\ . 
dy dx ) Sl 



sin 9 cos 9 



dUy 


du x \ 


, 1 | 


dx 


dy ) 


+ 2 < 


1 

+ 2 


( duy 


du x \ 


\ dy 


dx ) 



du,j du x 
dx + dy 

sin(20) . 



which 

(1.30) 



(1.31) 

(1.32) 

(1.33) 

(1.34) 

(1.35) 



cos(20) 



(1.36) 
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In the final expression (1.36), the first term is independent of the inclination 9 
of the line PQ and hence represents a rigid-body rotation as in Figure 1 .9(h). We 
denote this rotation by the symbol u, which with the convention illustrated in Figure 
1.8 is anticlockwise positive. 

Notice that as this term is independent of 9 in equation (1.36), it is the same for 
any right-handed set of axes — i.e. 



_ 1 (du y _ 1 / ddy ^ 

W 2 \ dx dy ) 2 y dx' dy' ) 1 



(1.37) 



for any x\y' . 

In three dimensions, u represents a small positive rotation about the 2 -axis and 
is therefore more properly denoted by ut 2 to distinguish it from the corresponding 
rotations about the x and y axes — i.e. 



_ 1 (du 2 du y \ _ 1 (du x du 2 \ 

Ux ~2 [dy ~ dz ) '' “ y ~2[dz dx) 



_lf duy _ du x \ 
Uz 2 \ dx dy ) 



(1.38) 



or in suffix notation 




(1.39) 



where is the alternating tensor which is defined to be 1 if the suffices are in cyclic 
order (e.g. 2,3,1), -1 if they are in reverse cyclic order (e.g. 2,1,3) and zero if any two 
suffices are the same. (Notice, the alternating tensor is not a second order Cartesian 
tensor.) 

The rotation a; is a vector in three-dimensional problems and can be defined by 
the equation 



1 , ^ 

u> = xcurl u = -V x u , 
2 2 



(1.40) 



but in two dimensions it behaves as a scalar since two of its components degenerate 
to zero 6 . 



1.2.4 Definition of shear strain 



We are now in a position to define the shear strain e xy as the difference between the 
rotation of a line drawn in the x-direction and the corresponding rigid-body rotation. 



utj, i.e. 



^xy 




1 / duy du x \ 

2 \ dx + dy ) 



(1.41) 



6 In some books, w zt w y ,uj z are denoted by tv yz ,uj z x,iv x ,j respectively. This notation is not used 
here because it gives the erroneous impression that u is a second order tensor rather than a vector. 
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and similarly 

1 (du z du y \ 1 (du x du z \ 

e ' JZ " 2 \dy dz) ' 6zx ~ 2 { dz dx ) ’ 



(1.42) 



Note that e xy so defined is one half of the quantity 7 xy used in Mechanics of 
Materials and in many older books on Elasticity. 

The strain-displacement relations (1.24, 1.25, 1.41, 1.42) can be written in the 
concise form 



dui duj 
dxj dx l 



(1.43) 



With the notation of equations (1.24, 1.25, 1.38, 1.41), we can now write (1.36) 
in the form 



<j)(PQ) =w 2 + e iy (cos 2 9 - sin 2 0) + (e^ - e xx ) sin 6 cos 0 0-44) 



and hence 

C x’y 1 — 



(by definition) 

= e iy (cos 2 9 - sin 2 #) + (e^ - e xx ) sin 0 cos 0 . (145) 

This is of course one of the coordinate transformation relations for strain. The 
other one 



^x’x’ — E xx cos 2 0 4- e yy sin 2 9 + 2e iy sin 9 cos 9 (1 .46) 

being obtainable from equations (1.24, 1.32, 1.33) in the same way. A comparison of 
equations (1.45, 1.46) and (1.7, 1.8) confirms that, with these definitions, the strain 
dj is a second order Cartesian tensor. 



1.3 Stress-strain relations 



We shall develop the various forms of the linear elastic stress-strain relations for the 
isotropic medium by regarding Young’s modulus E and Poisson’s ratio v as funda- 
mental constants. Hence we regard as experimentally determined the equations 



&XX 


&XX 

IT ' 


v °yy 
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1 'Ozz 
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(1-47) 
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a yy 
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1 'Ozz 

E 
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E 


(1.48) 


&ZZ — 


a zz 

E " 


vo xx 

E 


Wyy 

E ' 


(1.49) 



The relation between e xy and a xy can then be obtained by using the transformation 
relations. We know that there are three principal directions such that if we align them 
with x,y, z , we have 



xy ^yz — xx — 0 



(1.50) 
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and hence by symmetry 



*-xy — &yz — t-zx — 0 



(1.51) 



Using a coordinate system aligned with the principal directions, we write 

e*v = ( e yy - e xx ) sin 9 cos 9 (1.52) 



from equations (1.45, 1.51) 



(o yy - <T Iir )(l + i/) sin 9 cos 9 
E 

(1 + v)a X 'y' 

E 



from (1.47, 1.48, 1.8). 
We define 



A* = 



2(1 + t/) ’ 



so that equation (1.53) takes the form 



e *V 



_ a x‘y‘ 

2 jj. 



(1.53) 



(1.54) 



(1.55) 



1.3.1 Lame’s constants 

It is often desirable to solve equations (1.47-1.49) to express a xx in terms of e xx etc. 
The solution is routine and leads to the equation 



Ev(e xx 4 ~ Cyy + c xx ) Ee xx 
(i + v)[l-2v) (1 + ^) 



(1.56) 



and similar equations, which are more concisely written in the form 



etc., where 



and 



O xx — ^ C T 2 fJ>G X x 

Eu 2 fj,u 

(l + i/)(l-2 v) “ (1 - 2 1/) 



(1.57) 

(1.58) 



e = e xx + e yy + e zz = eu = div u (1.59) 

is known as the dilatation. 

The stress-strain equations (1.55, 1.57) can be written more concisely in the index 
notation in the form 



ffij — , 



(1.60) 
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where 5{j is the Kronecker delta, defined as 1 if i = j and 0 if i ^ j. The constants 
\,H are known as Lame’s constants. Young’s modulus and Poisson’s ratio can be 
written in terms of Lame’s constants through the equations 



/i(3A + 2/z) 
(A + n) 

A 

2(A + fi) 



(1.61) 

(1.62) 



1.3.2 Dilatation and bulk modulus 

The dilatation, e, is easily shown to be invariant as to coordinate transformation and 
is therefore a scalar quantity. In physical terms it is the local volumetric strain, since 
a unit cube increases under strain to a block of dimensions {l+e xx ), (l + e yy ), (l + e«) 
and hence the volume change is 

SV = (1 + e XI )(l + e yy )(l + e«) - 1 = e xx + e yy + e zz + 0{e xx e yy ) . (1.63) 

It can be shown that the dilatation e and the rotation vector u> are harmonic — 
i.e. V 2 e = V 2 u> = 0. For this reason, many early solutions of elasticity problems were 
formulated in terms of these variables, so as to make use of the wealth of mathematical 
knowledge about harmonic functions. We now have other more convenient ways of 
expressing elasticity problems in terms of harmonic functions, which will be discussed 
in Chapter 18 et seq.. 

The dilatation is proportional to the mean stress d through a constant known as 
the bulk modulus Kb- Thus, from equation (1.57), 

{&xx + Gyy + &zz) 1 ry /t za\ 

o = = ~(Tu = K b e , (1.64) 

where 

* i=A+ 5' ,= 3(r^r < L65 > 

We note that K b — > oo if u — > 0.5 — i.e. the material becomes incompressible. 



PROBLEMS 



1. Show that equations (1.47-1.49, 1.53) can be written in the concise form 



( 1 Y V^Oij VGmm&ij 



E 



E 



( 1 . 66 ) 



2. Prove that the partial derivatives d 2 f /dx 2 ', d 2 f/dxdy; d 2 f/dy' 2 of the scalar func- 
tion/ transform into the rotated coordinate system x\ y' by rules similar to equations 
(1.7-1. 9). 
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3. A rigid body is subjected to a small rotation u) z = fi 4^ 1 about the 2 -axis. If 
the displacement of the origin is zero, find expressions for the three displacement 
components u x , u v , u z asfunctions of x,y, z. 

4. Use the suffix notation to develop a general expression for the derivative 

dui 

dxj 



in terms of strains and rotations. 



5. Use the suffix notation to develop concise expressions for the three stress invariants 
h,h,h and the equivalent tensile stress ctr. 



6. Show that 



dxi 



(i) — — = Su and (ii) R = JxiXi , 
dxj 



where R — \R\ is the distance from the origin. Hence find dR/dxj in suffix notation. 
Confirm your result by finding dR/dx in x,y, z notation. 




Chapter 2 

EQUILIBRIUM AND 
COMPATIBILITY 



We can think of an elastic solid as a highly redundant framework — each particle is 
built-in to its neighbours. For such a framework, we expect to get some equations from 
considerations of equilibrium , but not as many as there are unknowns. The deficit 
is made up by compatibility conditions — statements that the deformed components 
must fit together. These latter conditions will relate the dimensions and hence the 
strains of the deformed components and in order to express them in terms of the same 
unknowns as the stresses (forces) we need to make use of the stress-strain relations 
as applied to each component separately. 

If we were to approximate the continuous elastic body by a system of intercon- 
nected elastic bars, this would be an exact description of the solution procedure. The 
only difference in treating the continuous medium is that the system of algebraic 
equations is replaced by partial differential equations describing the same physical or 
geometrical principles. 



2.1 Equilibrium equations 

We consider a small rectangular block of material — side 6x,Sy,Sz — as shown in 
Figure 2.1. We suppose that there is a body force 1 , p per unit volume and that 
the stresses vary with position so that the stress components on opposite faces of the 
block differ by the differential quantities Scr xx , 8o xy etc. Only those stress components 
which act in the x-direction are shown in Figure 2.1 for clarity. 



*A body force is one that acts directly on every particle of the body, rather than being applied 
by tractions at its boundaries and transmitted to the various particles by means of internal stresses. 
This is an important distinction which will be discussed further in Chapter 7, below. The commonest 
example of a body force is that due to gravity. 
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Figure 2.1: Forces in the x-direction on an elemental block. 

Resolving forces in the x-direction, we find 

(c^ii "b &O x x O xx }&y&Z -b ((7 xy ~b &0 X y 0 X y^SzSx 

+{o xz + Sa xt - a xz )5xSy + p x SxSySz = 0 . (2.1) 

Hence, dividing through by ( 5x5y5z ) and proceeding to the limit as these infinites- 
imals tend to zero, we obtain 



Similarly, we have 



do xx . do : xy 

dx dy 



dOyx d<7yy 
dx dy 



do x2 

+ -5 — + p x — 0 
dz 



+ ^ + P ,=0 



do zx do zy do 2Z 

dx dy dz 



Pz = 0 



( 2 . 2 ) 

(2.3) 

(2.4) 



or in suffix notation 



dcr a 

fej +Pi = 0 



These are the differential equations of equilibrium. 



2.2 Compatibility equations 

The easiest way to satisfy the equations of compatibility — as in framework problems 
— is to express all the strains in terms of the displacements. In a framework, this 
ensures that the components fit together by identifying the displacement of points in 
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two links which are pinned together by the same symbol. If the framework is redun- 
dant, the number of pin displacements thereby introduced is less than the number of 
component lengths (and hence extensions) determined by them — and the number 
of unknowns is therefore reduced. 

For a solid body, the process is essentially similar, but much more straightforward. 
We define the six components of strain in terms of displacements through equations 
(1.43). These six equations introduce only three unknowns (u x , u y , u z ) and hence the 
latter can be eliminated to give equations constraining the strain components. 

For example, from (1.24, 1.25) we find 



d 2 e xx _ d' i u x d 2 e yy _ d 3 u y 
dy 2 dxdy 2 ’ dx 2 dydx 2 

and hence 

d 2 e xx d 2 c yy _ d 2 t dux duj,\ _ d 2 e xy 

dy 2 dx 2 dxdy \dy dx ) dxdy ’ 

from (1.41) — i.e. 

d e xx _ e xy d e yy _ 

dy 2 dxdy dx 2 



( 2 . 6 ) 



(2.7) 



( 2 . 8 ) 



Two more equations of the same form may be obtained by permuting suffices. 
It is tempting to pursue an analogy with algebraic equations and argue that, since 
the six strain components are defined in terms of three independent displacement 
components, we must be able to develop three (i.e. 6-3) independent compatibility 
equations. However, it is easily verified that, in addition to the three equations similar 
to (2.8), the strains must satisfy three more equations of the form 



d 2 e zz _ d ( de yz de zx de xy \ 

dxdy dz\dx ^ dy dz ) 



(2.9) 



The resulting six equations are independent in the sense that no one of them can 
be derived from the other five, which all goes to show that arguments for algebraic 
equations do not always carry over to partial differential equations. 

A concise statement of the six compatibility equations can be written in the index 
notation in the form 



d ( foal _ _ n 

£pks dx/c \ dxi dxj J 



( 2 . 10 ) 



The full set of six equations makes the problem very complicated. In practice, 
therefore, most three-dimensional problems are treated in terms of displacements in- 
stead of strains. This satisfies the requirement of compatibility automatically. How- 
ever, in two dimensions, all except one of the compatibility equations degenerate to 
identities, so that a formulation in terms of stresses or strains is more practical. 
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2.2.1 The significance of the compatibility equations 

The physical meaning of equilibrium is fairly straightforward, but people often get 
mixed up about just what is being guaranteed by the compatibility equations. 

Single-valued displacements 

Mathematically, we might say that the strains are compatible when they are definable 
in terms of a single-valued, continuously differentiable displacement. 

We could imagine reversing this process — i.e. integrating the strains (displace- 
ment gradients) to find the relative displacement ( Ua~Ub ) of two points A, B in the 
solid (see Figure 2.2). 



A 

Figure 2.2: Path of the integral in equation (2.11). 
Formally we can write 2 




The integral will be along the line, S, and if the displacements are to be single-valued, 
it mustn’t make any difference if we change the line provided it remains within the 
solid. In other words, the integral should be path-independent. 

The compatibility equations are not quite sufficient to guarantee this. 

They do guarantee 

/ dS dS = ° ^ 2 ' 12 ^ 

around an infinitesimal closed loop. 

Now, we could make infinitesimal changes in our line from A to B by taking 
in such small loops until the whole line was sensibly changed, thus satisfying the 
requirement that the integral (2.11) is path-independent, but the fact that the line is 
changed infinitesimally stops us from taking a qualitatively (topologically) different 
route through a multiply-connected body. For example, in Figure 2.3, it is impossible 
to move Si to Sz by infinitesimal changes without passing outside the body 3 . 

Explicit forms of the integral (2.11) in terms of the strain components were developed by 
E.Cesaro and are known as Cesaro integrals. See, for example, A.E.H.Love, A Treatise on the 
Mathematical Theory of Elasticity, 4th edn.. Dover (1944), §156A. 

3 For a more rigorous discussion of this question, see A.E.H.Love, loc.cit. or B.A.Boley and 
J.H. Weiner, Theory of Thermal Stresses, John Wiley, New York, (1960), §§3.6-3. 8. 
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Figure 2.3: Qualitatively different integration paths in a multiply-connected body. 

In practice, if a solid is multiply-connected it is usually easier to work in terms of 
displacements and by-pass this problem. Otherwise the equivalence of topologically 
different paths in integrals like (2.11, 2.12) has to be explicitly enforced. 

Compatibility of deformed shapes 

A more ‘physical' way of thinking of compatibility is to state that the separate par- 
ticles of the body must deform under load in such a way that they fit together after 
deformation. This interpretation is conveniently explored by way of a ‘jig-saw’ anal- 
ogy- 

We consider a two-dimensional body cut up as a jig-saw puzzle (Figure 2.4), of 
which the pieces are deformable. Figure 2.4(a) shows the original puzzle and 2.4(6) 
the puzzle after deformation by some external loads Fj. (The pieces are shown as 
initially rectangular to aid visualization.) 

The deformation of the puzzle must satisfy the following conditions:- 

(i) The forces on any given piece (including external forces if any) must be in 
equilibrium. 

(ii) The deformed pieces must be the right shape to fit together to make the de- 
formed puzzle (Figure 2.4(6)). 

For the continuous solid, the compatibility condition guarantees that (ii) is satis- 
fied, since shape is defined by displacement derivatives and hence by strains. However, 
if the puzzle is multiply-connected — e.g. if it has a central hole — condition (ii) 
is not sufficient to ensure that the deformed pieces can be assembled into a coherent 
body. Suppose we imagine assembling the ‘deformed’ puzzle working from one piece 
outwards. The partially completed puzzle is simply-connected and the shape condi- 
tion is sufficient to ensure the success of our assembly until we reach a piece which 
would convert the partial puzzle to a multiply connected body. This piece may be 
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the right shape for both sides, but the wrong size for the separation. If so, it will be 
possible to leave the puzzle in a state with a discontinuity as shown in Figure 2.4(c) at 
any arbitrarily chosen position — i.e. there is no part of the body at which continuity 
breaks down — but there will be no way in which Figure 2A(b) can be constructed. 




Figure 2.4: A multiply-connected, deformable jig-saw puzzle: (a) Before deforma- 
tion; ( b ) After deformation; (c) Result of attempting to assemble the deformed puzzle 
if equation (2.12) is not satisfied for any closed path encircling the hole. 

The lines defining the two sides of the discontinuity in Figure 2.4(c) are the same 
shape and hence, in the most general case, the discontinuity can be defined by six 
arbitrary constants corresponding to the three rigid-body translations and three ro- 
tations needed to move one side to coincide with the other 4 . We therefore get six 
additional algebraic conditions for each hole in a multiply-connected body. 

4 This can be proved by evaluating the relative displacements of corresponding points on opposite 
sides of the cut, using an integral of the form (2.1 1) whose path does not cross the cut. 
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2.3 Equilibrium equations for displacements 



Remembering that it is often easier to work in terms of displacements to avoid com- 
plications with the compatibility conditions, it is convenient to express the eqilibrium 
equations in terms of displacements. This is done by substituting for the stresses 
from the stress-strain relations (1.60) into equation (2.2) giving 




+ 2/i 



/ de xx de xy de xz\ 

\ dx dy dz ) 



+ p x — 0 • 



(2.13) 



Also, 

2 f^i + ^y + <bEi\ 

\ ox dy dz J 



d 2 u x d 2 u x d 2 u y d 2 u x d 2 u z 
dx 2 + dy 2 + dxdy + dz 2 + dxdz 

V2 “ i + £’ (2 ' 14) 



and hence equation (2.13) can be written 



(A + /i) — div u + pV 2 Ux + Pi = 0 . 



(2.15) 



Two similar equations are obtained for the equilibrium conditions in the y- and 
^-directions. 

We also note that 

2v 



X 



+ /i = /i ( 1 + rrk)-(T^ 



(1 - 2 «,) ’ 



(2.16) 



from (1.58) and hence the general equilibrium condition can be concisely written in 
the vector form 

Vdiv u + (1 - 2u)V 2 u + ^ ~ 2 ^ P - = 0 , (2.17) 

or in suffix notation 



+ ( i-2„) 

dXidXj V j dx k dx k P 



= 0 . 



(2.18) 



PROBLEMS 

1. Show that, if there are no body forces, the dilatation e must satisfy the condition 

V 2 e = 0 . 

2. Show that, if there are no body forces, the rotation u> must satisfy the condition 

V 2 u> = 0 . 
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3. One way of satisfying the compatibility equations in the absence of rotation is to 
define the components of displacement in terms of a potential function ip through the 
relations 





u z 



dip 

dz 



Use the stress-strain relations to derive expressions for the stress components in 
terms of ip. 

Hence show that the stresses will satisfy the equilibrium equations in the absence 
of body forces if and only if 

V 2 ip = constant . 



4. Plastic deformation during a manufacturing process generates a state of residual 
stress in the large body z > 0. If the residual stresses are functions of z only and the 
surface z = 0 is not loaded, show that the stress components a vz , a zx , a zz must be zero 
everywhere. 

5. By considering the equilibrium of a small element of material similar to that shown 
in Figure 1.2, derive the three equations of equilibrium in cylindrical polar coordinates 
r,6,z. 

6. In cylindrical polar coordinates, the strain-displacement relations for the ‘in-plane’ 
strains are 

du T 1 / 1 du r dug u e\ u r , 1 dug 

Crr dr ' ^ 2 86 + dr r ) ' Cae r rdO 

Use these relations to obtain a compatibility equation that must be satisfied by the 
three strains. 



7. If no stresses occur in a body, an increase in temperature T causes unrestrained 
thermal expansion defined by the strains 



&XZ Cyy — 6 ZZ — CtT | exy 6y Z — € zx — 0 . 

Show that this is possible only if T is a linear function of x, y, z and that otherwise 
stresses must be induced in the body, regardless of the boundary conditions. 



8. If there are no body forces, show that the equations of equilibrium and compati- 
bility imply that 



( 1 +") 



9 * 01 , 



+ 



d 2 oi 



kk 



dxkdx k dx{dxj 



= 0 
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TWO-DIMENSIONAL 

PROBLEMS 




Chapter 3 

PLANE STRAIN AND PLANE 
STRESS 



A problem is two-dimensional if the field quantities such as stress and displacement 
depend on only two coordinates ( x,y ) and the boundary conditions are imposed on 
a line f(x,y) = 0 in the xy-plane. 

In this sense, there are strictly no two-dimensional problems in elasticity. There 
are circumstances in which the stresses are independent of the z-codrdinate, but all 
real bodies must have some bounding surfaces which are not represented by a line 
in the xy-plane. The two-dimensionality of the resulting fields depends upon the 
boundary conditions on such surfaces being of an appropriate form. 



3.1 Plane strain 

It might be argued that a closed line in the xy-plane does define a solid body — namely 
an infinite cylinder of the appropriate cross-section whose axis is parallel to the z- 
direction. However, making a body infinite does not really dispose of the question of 
boundary conditions, since there are usually some implied boundary conditions ‘at 
infinity’. For example, the infinite cylinder could be in a state of uniaxial tension, 
er 22 = C, where C is an arbitrary constant. However, a unique two-dimensional infinite 
cylinder problem can be defined by demanding that u x , u y be independent of z and 
that it 2 =0 for allx,y, z, in which case it follows that 

^ zx — ^zy ^zz — 0 • ( 3 - 1 ) 

This is the two-dimensional state known as plane strain. 

In view of the stress-strain relations, an equivalent statement to equation (3.1) is 

Ozx — a zy = 0 ! u z = 0 , (3.2) 

and hence a condition of plane strain will exist in a finite cylinder provided that (i) 
any tractions or displacements imposed on the sides of the cylinder are independent 
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of z and have no component in the z-direction and (ii) the cylinder has plane ends 
(e.g. z = ±c) which are in frictionless contact with two plane rigid walls. 

From the condition e zz =0 (3.1) and the stress-strain relations, we can deduce 



— i.e. 
and hence 



O zi v(Oxx 3 dyy ) 

~E E 



°zz — V\Oxx 


+ Oyy) 


O XT 


v(J yy 




IT ~ 


E 


E 


(1 - u 2 )a xx 


u(l + V 


E 


E 



(3-3) 



(3.4) 



Of course, there are comparatively few practical applications in which a cylinder 
with plane ends is constrained between frictionless rigid walls, but fortunately the 
plane strain solution can be used in an approximate sense for a cylinder with any 
end conditions, provided that the length of the cylinder is large compared with its 
cross-sectional dimensions. 

We shall illustrate this with reference to the long cylinder of Figure 3.1, for which 
the ends, z = ±c are traction-free and the sides are loaded by tractions which are 
independent of z. 




We first solve the problem under the plane strain assumption, obtaining an exact 
solution in which all the stresses are independent of z and in which there exists a 
normal stress a zz on all z-planes, which we can calculate from equation (3.3). 
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The plane strain solution satisfies all the boundary conditions of the problem 
except that a zz also acts on the end faces, z — ±c, where it appears as an unwanted 
normal traction. We therefore seek a corrective solution which, when superposed on 
the plane strain solution, removes the unwanted normal tractions on the surfaces 
z = ±c, without changing the boundary conditions on the sides of the cylinder. 

3.1.1 The corrective solution 

The corrective solution must have zero tractions on the sides of the cylinder and a 
prescribed normal traction (equal and opposite to that obtained in the plane strain 
solution) on the end faces. 

This is a fully three-dimensional problem which generally has no closed-form so- 
lution. However, if the prescribed tractions have the linear form 

o zz = B + Cx + Dy , ( 3 . 5 ) 

the solution can be obtained in the context of Mechanics of Materials by treating the 
long cylinder as a beam subjected to an axial force 

F = J A Cr 2Z dA ( 3 . 6 ) 

and bending moments 

M x = a zz ydA ; M y — - a zz xda ( 3 . 7 ) 

J A J A 

about the axes Ox, Oy respectively, where 0 is chosen to coincide with the centroid 
of the cross-section 1 . 

If the original plane strain solution does not give a distribution of o zz of this 
convenient linear form, we can still use equations (3.6, 3.7) to define the force and 
moments for an approximate Mechanics of Materials corrective solution. The error 
involved in using this approximate solution will be that associated with yet another 
corrective solution corresponding to the problem in which the end faces of the cylinder 
are loaded by tractions equal to the difference between those in the plane strain 
solution and the Mechanics of Materials linear form (3.5) associated with the force 
resultants (3.6, 3.7). 

In this final corrective solution, the ends of the cylinder are loaded by self- 
equilibrated tractions , since the Mechanics of Materials approximation is carefully 
chosen to have the same force and moment resultants as the required exact solution 
and in such cases, we postulate that significant stresses will only be generated in the 
immediate vicinity of the ends — or more precisely, in regions whose distance from the 

*The Mechanics of Materials solution for axial force and pure bending is in fact exact in the sense 
of the Theory of Elasticity, but we shall not prove this here. See for example S.P. Timoshenko and 
J.N.Goodier, Theory of Elasticity, McGraw-Hill, New York, 3rd. edn. (1970) §102. 
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ends is comparable with the cross-sectional dimensions of the cylinder. If the cylinder 
is many times longer than its cross-sectional dimensions, there will be a substantial 
portion near the centre where the final corrective solution gives negligible stresses 
and hence where the sum of the original plane strain solution and the Mechanics of 
Materials correction is a good approximation to the actual three-dimensional stress 
field. 



3.1.2 Saint-Venant’s principle 

The thesis that a self-equilibrating system of loads produces only local effects is known 
as Saint-Venant’s principle. It seems intuitively reasonable, but has not been proved 
rigorously except for certain special cases — some of which we shall encounter later 
in this book (see for example Chapter 6). It can be seen as a consequence of the rule 
that alternate load paths through a structure share the load in proportion with their 
stiffnesses. If a region of the boundary is loaded by a self-equilibrating system of 
tractions, the stiffest paths are the shortest — i.e. those which do not penetrate far 
from the loaded region. Hence, the longer paths — which are those which contribute 
to stresses distant from the loaded region — carry relatively little load. 

Note that if the local tractions are not self-equilibrating, some of the load paths 
must go to other distant parts of the boundary and hence there will be significant 
stresses in intermediate regions. For this reason, it is important to superpose the 
Mechanics of Materials approximate corrective solution when solving plane strain 
problems for long cylinders with traction free ends. By contrast, the final stage of 
solving the ‘Saint-Venant problem to calculate the correct stresses near the ends is 
seldom of much importance in practical problems, since the ends, being traction-free, 
are not generally points of such high stress as the interior. 

As a point of terminology, we shall refer to problems in which the boundary condi- 
tions on the ends have been corrected only in the sense of force and moment resultants 
as being solved in the weak sense with respect to these boundaries. Boundary condi- 
tions are satisfied in the strong sense when the tractions are specified in a pointwise 
rather than a force resultant sense. 



3.2 Plane stress 

Plane stress is an approximate solution, in contrast to plane strain, which is exact. 
In other words, plane strain is a special solution of the complete three-dimensional 
equations of elasticity, whereas plane stress is only approached in the limit as the 
thickness of the loaded body tends to zero. 

It is argued that if the two bounding 2 -planes of a thin plate are sufficiently close 
in comparison with the other dimensions, and if they are also free of tractions, the 
stresses on all parallel 2-planes will be sufficiently small to be neglected in which case 
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we write 

& zx — G zy ~~ 0 zz = 0 » 

for all x,y , z. 

It then follows that 



but e«=^0, being given in fact by 

v . . 

e zz — O \&xx T a yy) • 

The two-dimensional stress-strain relations are then 

Gxx Wyy 

~E E~ 

Oyy V&xx 



&XX — 

e yy = 



( 3 . 8 ) 

( 3 . 9 ) 

( 3 . 10 ) 

( 3 . 11 ) 

( 3 . 12 ) 



The fact that plane stress is not an exact solution can best be explained by con- 
sidering the compatibility equation 



d 7 e yy d 2 e yz d 2 e zz 

dz 2 dydz dy 2 



( 3 . 13 ) 



Since ex-hypothesi none of the stresses vary with z, the first two terms in this 
equation are identically zero and hence the equation will be satisfied if and only if 



d 2 e zz 

dy 2 



= 0 . 



( 3 . 14 ) 



This in turn requires 

d ip XX ~b <7yy) _ p. 

dy 2 

from equation (3.10). 

Applying similar arguments to the other compatibility equations, we conclude 
that the plane stress assumption is exact if and only if [p xx +a yy ) is a linear function 
of X,y. i.e. if 

G XX T tTyy — 13 ~ 1“ Cx + Dy. (3.16) 

The attentive reader will notice that this condition is exactly equivalent to (3.3, 
3.5). In other words, the plane stress solution is exact if and only if the ‘weak form’ 
solution of the corresponding plane strain solution is exact. Of course this is not 
a coincidence. In this special case, the process of off-loading the ends of the long 
cylinder in the plane strain solution has the effect of making a zz zero throughout the 
cylinder and hence the resulting solution also satisfies the plane stress assumption, 
whilst remaining exact. 
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3.2.1 Generalized plane stress 



The approximate nature of the plane stress formulation is distasteful to elasticians 
of a more mathematical temperament, who prefer to preserve the rigour of an exact 
theory. This can be done by the contrivance of defining the average stresses across 
the thickness of the plate — e.g. 



*-kf. 



& xxdcC * 



(3.17) 



It can then be shown that the average stresses so defined satisfy the plane stress 
equations exactly. This is referred to as the generalized plane stress formulation. 

In practice, of course, the gain in rigour is illusory unless we can also establish that 
the stress variation across the section is small, so that the local values are reasonably 
close to the average. A fully three-dimensional theory of thin plates under in-plane 
loading shows that the plane stress assumption is a good approximation except in 
regions whose distance from the boundary is comparable with the plate thickness 2 . 



3.2.2 Relationship between plane stress and plane strain 

The solution of a problem under either the plane strain or plane stress assumptions 
involves finding a two-dimensional stress field, defined in terms of the components 
o xx , &xyi °yy> which satisfies the equilibrium equations (2.5), and for which the corre- 
sponding strains, e XXl e xy , e yy , satisfy the only non-trivial compatibility equation (2.8). 
The equilibrium and compatibility equations are the same in both formulations, the 
only difference being in the relation between the stress and strain components, which 
for normal stresses are given by (3.4) for plane strain and (3.11, 3.12) for plane stress. 
The relation between the shear stress o xy and the shear strain e xy is the same for both 
formulations and is given by equation (1.55). Thus, from a mathematical perspective, 
the plane strain solution simply looks like the plane stress solution for a material with 
different elastic constants and vice versa. In fact, it is easily verified that equation 
(3.4) can be obtained from (3.11) by making the substitutions 

F= & S 

(1-i/ 2 ) ’ " (l-i/) 



and then dropping the primes. This substitution also leaves the shear stress-shear 
strain relation unchanged as required. 

An alternative approach is to write the two-dimensional stress-strain relations in 
the form 







/c — f- 1 
8/i 



&XX 



3 — K 

8/i 



'yy 



k + 1 

8/i 



8/i 



e iy 



_ u xy 



2/i 

(3.19) 



2 See for example S.P.Timoshenko and J.N.Goodier, toe. cit., §98. 
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where k — 3-4^ for plane strain and (3 - u)/{\ + u) for plane stress and is known as 
Kolosov’s constant. 

In the following chapters, we shall generally treat two-dimensional problems with 
the plane stress assumptions, noting that results for plane strain can be recovered by 
the substitution (3.18) when required. 

PROBLEMS 



1. The plane strain solution for the stresses in the rectangular block 0 <x<a,-b< 
y<b,—c<z<c with a given loading is 

3 Fxy 3 .F(6 2 - y 2 ) _ 3 vFxy 

° xx ~ ~2&~ ’ ° xy ~ 463 ’ Cyy ~ 0 ' ° lz ~ 2p~~ - 

Find the tractions on the surfaces of the block and illustrate the results on a sketch 
of the block. 

We wish to use this solution to solve the corresponding problem in which the 
surfaces z = ic are traction-free. Determine an approximate corrective solution for 
this problem by offloading the unwanted force and moment resultants using the ele- 
mentary bending theory. Find the maximum error in the stress a zz in the corrected 
solution and compare it with the maximum tensile stress in the plane strain solution. 



2. For a solid in a state of plane stress, show that if there are body forces p x ,p y per 
unit volume in the direction of the axes x, y respectively, the compatibility equation 
can be expressed in the form 



V 2 (a xx +a yy ) = -(1 + v) 



(dVx dp y \ 

{ dx + dy ) • 



Hence deduce that the stress distribution for any particular case is independent 
of the material constants and the body forces, provided the latter are constant. 



3. (i) Show that the compatibility equation (2.8) is satisfied by unrestrained thermal 
expansion ( e xx = e yy = aT,e xy = 0), provided that the temperature, T, is a 
two-dimensional harmonic function — i.e. 

d 2 T d 2 T 
dx 2 + dy 2 

(ii) Hence deduce that, subject to certain restrictions which you should explicitly 
list, no thermal stresses will be induced in a thin body with a steady-state, 
two-dimensional temperature distribution and no boundary tractions. 

(iii) Show that an initially straight line on such a body will be distorted by the heat 
flow in such a way that its curvature is proportional to the local heat flux across 
it. 
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4. Find the inverse relations to equations (3.18) — i.e. the substitutions that should 
be made for the elastic constants E, u in a plane strain solution if we want to recover 
the solution of the corresponding plane stress problem. 



5. Show that in a state of plane stress without body forces, the in-plane displacements 
must satisfy the equations 



V 2 u x + 



(1 + v) d 
(1 - v) dx 



V dx dy ) 



= 0 



V z u y + 



(l + i/) d 
(1 - u) dy 



fdux du y \ 
V dx dy) 



= 0 . 



6. Show that in a state of plane strain without body forces, 

de _ /I — 2//N du z de __ /I - 2tA du> z 
dx \ 1 - i ' ) dy ' dy \ 1 — v ) dx 



7. If a material is incompressible (i/ = 0.5), a state of hydrostatic stress o xx —o yy — o zz 
produces no strain. One way to write the corresponding stress-strain relations is 

Oij — 2/rey — (]S{j , 

where q is an unknown hydrostatic pressure which will generally vary with position. 
Also, the condition of incompressibility requires that the dilatation 

e = e kk = 0 . 

Show that the stress components and the hydrostatic pressure q must satisfy the 
equations 

V 2 q = div p and a xx + a yy = -2q , 
where p is the body force. 




Chapter 4 



STRESS FUNCTION 
FORMULATION 

4.1 The concept of a scalar stress function 



Newton’s law of gravitation states that two heavy bodies attract each other with a 
force proportional to the inverse square of their distance — thus it is essentially a 
vector theory, being concerned with forces. However, the idea of a scalar gravitational 
potential can be introduced by defining the work done in moving a unit mass from 
infinity to a given point in the field. The principle of conservation of energy requires 
that this be a unique function of position and it is easy to show that the gravitational 
force at any point is then proportional to the gradient of this scalar potential. Thus, 
the original vector problem is reduced to a problem about a scalar potential and its 
derivatives. 

In general, scalars are much easier to deal with than vectors. In particular, they 
lend themselves very easily to coordinate transformations, whereas vectors (and to an 
even greater extent tensors) require a set of special transformation rules (e.g. Mohr’s 
circle). 

In certain field theories, the scalar potential has an obvious physical significance. 
For example, in the conduction of heat, the temperature is a scalar potential in terms 
of which the vector heat flux can be defined. However, it is not necessary to the 
method that such a physical interpretation can be given. The gravitational potential 
can be given a physical interpretation as discussed above, but this interpretation may 
never feature in the solution of a particular problem, which is simply an excercise 
in the solution of a certain partial differential equation with appropriate boundary 
conditions. In the theory of elasticity, we make use of scalar potentials called stress 
functions or displacement functions which have no obvious physical meaning other 
than their use in defining stress or displacement components in terms of derivatives. 
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4.2 Choice of a suitable form 

In the choice of a suitable form for a stress or displacement function, there is only one 
absolute rule — that the operators which define the relationship between the scalar 
and vector (or tensor) quantities should indeed define a vector (or tensor). 

For example, it is appropriate to define the displacement in terms of the first 
derivatives (the gradient) of a scalar or to define the stress components in terms of 
the second derivatives of a scalar, since the second derivatives of a scalar form the 
components of a Cartesian tensor. 

In effect, what we are doing in requiring this similarity of form between the def- 
initions and the defined quantity is ensuring that the relationship is preserved in 
coordinate transformations. It would be quite possible to work out an elasticity prob- 
lem in terms of the displacement components u x ,u y ,u z , treating these as essentially 
scalar quantities which vary with position — indeed this was a technique which was 
used in early theories. However, we would then get into trouble as soon as we tried 
to make any statements about quantities in other coordinate directions. By contrast, 
if we define (for example) u x = dip/dx\ u y = dfjdy, u z — dtp/dz, (i.c. u - V(/i) , it 
immediately follows that u' x = dip/dx' for any x' . 



4.3 The Airy stress function 



The two-dimensional problem of elasticity is most easily reduced to a tractable po- 
tential problem by representing the stress components in the form 

d 2 (j> d 2 <t> _ d 2 <p 

W ; ° w = fa 2 ' ° xy ~ ' 



&XX — 



(4.1) 



This representation was introduced by G.B.Airy 1 in 1862 and <p is therefore gen- 
erally referred to as the Airy stress function. It is not the most obvious form. It 
would seem more natural to write a xx — d 2 ip/dx 2 \ o yy = d 2 tp/dy 2 \ o xy = d 2 i/>/dxdy 
and indeed this also leads to a representation which is widely used as part of the 
general three-dimensional solution (see Chapter 18 below). 

It is easily verified that equation (4.1) transforms as a Cartesian tensor as required. 
For example, using (1.32) we can write 



d 2 <p 

dx’ 2 



„ d d 

cos 9— + sinfl— 
ox oy 



<t> 



8 2 (j) 2 „ 9 2 l i> .o n n d 2< t> . 

—fr cos 2 6 + —fr sin 2 6 + 2— 4- sin 6 cos 0 
ox i ay 1 oxoy 



(4.2) 



from which using (1.9 and 4.1) we deduce that o,y y ' = d 2 (t)/ dx 12 as required. 

'For a good historical survey of the development of potential function methods in Elasticity, see 
H.M.Westergaard, Theory of Elasticity and Plasticity, Dover, New York (1964), Chapter 2. 
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4.4 The governing equation 

As discussed in Chapter 2, the stress or displacement field has to satisfy the equations 
of equilibrium and compatibility if they are to describe permissible states of an elastic 
body. In stress function representations, this generally imposes certain constraints on 
the choice of stress function, which can be expressed by requiring it to be a solution of 
a certain partial differential equation. We shall determine this governing equation by 
substituting the representation (4.1) into the equilibrium and compatibility equations 
in two dimensions. 



4.4.1 The equilibrium equations 

On substituting (4.1) into the first equilibrium equation, we obtain 

d 3 <}> _ d 3 <t> 
dy 2 dx dxdy 2 ^ 



(4.3) 



and hence p x = 0. In the same way, the second equilibrium equation is satisfied as 
long as p y = 0. Thus, the Airy stress function automatically satisfies the equilibrium 
equations provided the body forces are zero. In some books, it is stated that this is 
precisely the reason for preferring the Airy function to the more obvious form referred 
to in §4.3 above, but this is misleading, since the ‘more obvious form’ — although 
it requires some constraints on ip in order to satisfy equilibrium — can be shown to 
define displacements which automatically satisfy the compatibility condition 2 . The 
real reason for preferring the Airy function is that it is capable of describing all 
possible two-dimensional states of stress in an elastic body, whilst the alternative is 
more restrictive. This will not be proved here. 

It is worth noting that the Airy stress function can be used in other applica- 
tions, with inelastic constitutive laws — for example plasticity theory — where its 
satisfaction of the equilibrium equations remains an advantage. 



4.4.2 Non-zero body forces 

If the body force p is not zero, but is of a restricted form such that it can be written 
p = — W, where V is a scalar potential, equations (4.1) can be generalized by 
including an extra term +V in each of the normal stress components whilst leaving 
the shear stress definition unchanged. It is easily verified that, with this modification, 
the equilibrium equations are again satisfied. 

Problems involving body forces will be discussed in more detail in Chapter 7 
below. For the moment, we restrict attention to the the case where p — 0 and hence 
where the representation (4.1) is appropriate. 

2 See Problem 2.3. 
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4.4.3 The compatibility condition 



In order to determine the constraints placed upon the choice of 4> by the compati- 
bility condition, we first express the latter in terms of stresses using the stress-strain 
relations, obtaining 



d 2 a. 



2 ^2 _ 

” - - 2(1 + vy-^L 



+ 



d 2 o. 



vv _ L P ffxx _ Q 



dy 2 dy 2 dxdy dx 2 dx 2 

We then substitute for the stress components from (4.1) obtaining 

^ d*<}> 



ay _ 

c?y 4 V dx 2 dy 2 



d i 6 , d^cb 

V + 2(1 + i/) V 



dx 2 dy 2 + ax 4 V dx 2 dy 2 



= 0 



i.e. 



ay 2 av 

ax 4 dx 2 dy 2 



ay . _ _ 

+ ^ " 1 ^72 + ^2 j ^ -0 



ax 2 ay 2 



(4.4) 



(4.5) 



(4.6) 



This equation is known as the biharmonic equation and is usually written in the 
concise form 

vy = 0 . 



The biharmonic equation is the governing equation for the Airy stress function 
in elasticity problems. Thus, by using the Airy stress function representation, the 
problem of determining the stresses in an elastic body is reduced to that of finding 
a solution of equation (4.7) (i.e. a biharmonic function ) whose derivatives satisfy 
certain boundary conditions on the surfaces. 



4.4.4 Method of solution 

Historically, the boundary-value problem for the Airy stress function has been ap- 
proached in a semi-inverse way — i.e. by using the variation of tractions along the 
boundaries to give a clue to the kind of function required, but then exploring the 
stress fields developed from a wide range of such functions and selecting a combina- 
tion which can be made to satisfy the required conditions. The disadvantage with 
this method is that it requires a wide experience of particular solutions and even then 
is not guaranteed to be successful. 

A more modern method which has the advantage of always developing an appro- 
priate stress function if the boundary conditions are unmixed (i.e. all specified in 
terms of stresses or all in terms of displacements) is based on representing the stress 
function in terms of analytic functions of the complex variable. This method will 
be discussed briefly in Chapter 19. However, although it is powerful and extremely 
elegant, it is an unfortunate fact that for most problems it leaves us with the task of 
evaluating a difficult contour integral, which may not be as convenient a numerical 
technique as a direct series or finite difference attack on the original problem using 
real stress functions. (There are some exceptions — notably for bodies which are 
susceptible of a simple conformal transformation into the unit circle.). 
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4.4.5 Reduced dependence on elastic constants 

It is clear from dimensional considerations that, when the compatibility equation is 
expressed in terms of 0, Young’s modulus must appear in every term and can therefore 
be cancelled, but it is an unexpected bonus that the Poisson’s ratio terms also cancel 
in (4.5), leaving an equation that is independent of elastic constants. It follows 
that the stress field in a simply-connected elastic body 3 in a state of plane strain 
or plane stress is independent of the material properties if the boundary conditions 
are expressed in terms of tractions and in particular, that the plane stress and plane 
strain fields are identical. 

Dundurs 4 has shown that a similar reduced dependence on elastic constants occurs 
in plane problems involving interfaces between two dissimilar elastic materials. In 
such cases, three independent dimensionless parameters (e.g. M 1 /M 2 , ^l, ^ 2 ) can be 
formed from the elastic constants Mi, ^l, M 2 , ^2 of the materials 1, 2 respectively, but 
Dundurs proved that the stress field can be written in terms of only two parameters 
which he defined as 

_ / Ki + 1 K2 + l^ //«! + 1 «2"bl\ 

\ Mi M 2 J / \ Mi M 2 / 

P - ( Kl ~ 1 *2 ~ 1 \ /V «1 + 1 | «2 + 1 \ 

\ Ml M 2 / / \ Ml M 2 / ’ 

where k, is Kolosov's constant (see equation (3.19)). It can be shown that Dundurs' 
parameters must lie in the range — l<a< 1 ; — 0.5 < /? < 0.5 if 0 <V\, V 2 <0.5. 



(4.8) 

(4.9) 



PROBLEMS 



1. Newton's law of gravitation states that two heavy particles of mass mi, m 2 respec- 
tively will experience a mutual attractive force 

_ 7^i m 2 

F? ' 



where R is the distance between the particles and 7 is the gravitational constant. Use 
an energy argument and superposition to show that the force acting on a particle of 
mass mo can be written 

F = ~7mo VU , 



where 






m(£, 7, O d tdr/d( 



O 2 + (y - v ) 2 + (* - O 2 



3 The restriction to simply-connected bodies is necessary, since in a multiply-connected body there 
is an implied displacement condition, as explained in §2.2 above. 

4 J. Dundurs. Discussion on ‘Edge bonded dissimilar orthogonal elastic wedges under normal and 
shear loading, ASME J.Appl.Mech., Vol. 36 (1969), 650-652. 
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represents the volume of the universe and p is the density of material in the universe, 
which will generally be a function of position (£, T), Q. 

Could a similar method have been used if Newton's law had been of the more 
general form 

7 m 1 m 2 
F ‘ 

If so, what would have been the corresponding expression for VI If not, why not? 

2. An ionized liquid in an electric field experiences a body force p. Show that the 
liquid can be in equilibrium only if p is a conservative vector field. Hint: Remember 
that a stationary liquid must be everywhere in a state of hydrostatic stress. 

3. (i) Show that the function 

<t> = yu) + ip 

satisfies the biharmonic equation provided that ui,ip are both harmonic (i.e. 

V 2 w = 0,V 2 tf< = 0). 



(ii) Develop expressions for the stress components in terms of u>, ip, based on the 
use of ip as an Airy stress function. 



(iii) Show that a solution suitable for the half-plane y>0 subject to normal surface 
tractions only (i.e. <7 Iy = 0 on y = 0)can be obtained by writing 

dip 




and hence that under these conditions the normal stress o xx near the surface 
y = 0 is equal to the applied traction a yy . 

(iv) Do you think this is a rigorous proof? Can you think of any exceptions? If so, 
at what point in your proof of section (iii) can you find a lack of generality? 



4. The constitutive law for an orthotropic elastic material in plane stress can be 
written 



C X x — ^ll^xx "b ^12 &yy i ^ yy — S\2&xx T $22&yy I &xy — S44G xy , 

where Sn,Si 2 ,S 22 -, S 44 are elastic constants. 

Using the Airy stress function <p to represent the stress components, find the 
equation that must be satisfied by <p. 

5. Show that if the two-dimensional function ui{x,y) is harmonic (V 2 w = 0), the 
function 

(p-{x 2 + y 2 )u) 
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will be biharmonic. 

6. The constitutive law for an incompressible elastic material can be written 

@ij tx &ij T 2 /r e j j , 



where 



a — 



Okk 

3 



represents an arbitrary hydrostatic stress field. Some soils can be approximated as 
an incompressible material whose modulus varies linearly with depth, so that 



/r = Mz 



for the half space z> 0. 

Use the displacement function representation 



u — V(f> 



to develop a potential function solution for the stresses in such a body. Show that 
the functions <j>, a must satisfy the relations 

V'V = 0 ; o = -2 M d 4- 
az 

and hence obtain expressions for the stress components in terms of the single harmonic 
function cj). 

If the half-space is loaded by a normal pressure 

cr Z z{x,y,0) = -p{x y y) ; a zx {x,y,0) = a zy (x,y,Q) = 0 , 

show that the corresponding normal surface displacement u z {x,y, 0) is linearly pro- 
portional to the local pressure p{x, y) and find the constant of proportionality 5 . 

7. Show that Dundurs’ constant /?-> 0 for plane strain in the limit where =0.5 and 
y.\/p 2 — >0 — i.e. material 1 is incompressible and has a much lower shear modulus 6 
than material 2. What is the value of a in this limit? 

8. Solve Problem 3.7 for the case where there is no body force, using the Airy 
stress function 0 to represent the stress components. Hence show that the governing 
equation is V‘ 1 i/> = 0, as in the case of compressible materials. 

5 For an alternative proof of this result see C.R.Calladine and J. A. Greenwood, Line and point 
loads on a non-homogeneous incompressible elastic half-space, Q.J.Mech.Appl.Math., Vol 31 (1978), 
507-529. 

s This is a reasonable approximation for the important case of rubber (material 1) bonded to steel 
(material 2). 




Chapter 5 

PROBLEMS IN RECTANGULAR 
COORDINATES 



The Cartesian coordinate system (a:, y ) is clearly particularly suited to the problem 
of determining the stresses in a rectangular body whose boundaries are defined by 
equations of the form x = a,y = b. A wide range of such problems can be treated using 
stress functions which are polynomials in x, y. In particular, polynomial solutions can 
be obtained for ‘Mechanics of Materials’ type beam problems in which a rectangular 
bar is bent by an end load or by a distributed load on one or both faces. 



5.1 Biharmonic polynomial functions 

In rectangular coordinates, the biharmonic equation takes the form 

(5. 

and it follows that any polynomial in x, y of degree less than four will be biharmonic 
and is therefore appropriate as a stress function. However, for higher order polynomial 
terms, equation (5.1) is not identically satisfied. Suppose, for example, that we 
consider just those terms in a general polynomial whose combined degree (the sum 
of the powers of x and y) is N. We can write these terms in the form 

Pn{x, y) = A 0 x N + A x x N ~ l y + A 2 x N ~ 2 y 2 + . . . + A N y N (5.2) 

= (5.3) 

i=0 

where we note that there are {N + 1) independent coefficients, <4*(z = 0, N). If we 
now substitute P^(x,y) into equation (5.1), we shall obtain a new polynomial of 
degree (N - 4), since each term is differentiated four times. We can denote this new 



ay ay ay 

dx 4 + dx 2 dy 2 + dy 4 
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polynomial by Qn- 4 (x,y) where 

Qn-\{x, y) = S7 4 P N {x,y) (5.4) 

= £ . (5.5) 

t— o 

The (IV- 3) coefficients B 0 , . . . ,Bs-\ are easily obtained by expanding the right 
hand side of equation (5.4) and equating coefficients. For example, 

B 0 = N(N - l)(N - 2){N - 3Mo + 4(IV - 2){N - 3 )A 2 + 24A< . (5.6) 

Now the original function P^(x,y) will be biharmonic if and only if Qyv- 4 (z,t/) 
is zero for all x,y and this in turn is only possible if every term in the series (5.5) is 
identically zero, since the polynomial terms are all linearly independent of each other. 
In other words 

Bi = 0 ; i = 0 to W - 4 . (5.7) 

These conditions can be converted into a corresponding set of (N - 3) equations 
for the coefficients A{. For example, the equation jB o = 0 gives 

N(N - 1)(N - 2 )(N - 3)A> + 4 {N - 2 )(N - 3 )A 2 + 24A 4 = 0 , (5.8) 

from equation (5.6). We shall refer to the (N - 3) equations of this form as constraints 
on the coefficients Ai, since the coefficients are constrained to satisfy them if the 
original polynomial is to be biharmonic. 

One approach would be to use the constraint equations to eliminate (N - 3) of 
the unknown coefficients in the original polynomial — for example, we could treat 
the first fourcoefficients, Aq, Ai, A%, A 3 , as unknown constants and use the constraint 
equations to define all the remaining coefficients in terms of these unknowns. Equation 
(5.8) would then be treated as an equation for A 4 and the subsequent constraint 
equations would each define one new constant in the series. It may help to consider 
a particular example at this stage. Suppose we consider the fifth degree polynomial 



y) = A)Z 5 + A\X 4 y + A 2 x 3 y 2 + A 3 x 2 y 3 + A 4 xy 4 + A 5 y 5 , (5.9) 

which has six independent coefficients. Substituting into equation (5.4), we obtain 
the first degree polynomial 

Qi(x,y) = (12(L4 0 + 244 2 + 24 A 4 )x + (24 Ai + 24 A 3 + l20A 5 )y . (5.10) 

The coefficients of x and y in must both be zero if P 5 is to be biharmonic and 
we can write the resulting two constraint equations in the form 

A a = -5A 0 -A 2 (5.11) 

A 5 = -A,/5-A 3 /5. (5.12) 
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Finally, we use (5.11, 5.12) to eliminate A 4 ,A 5 in the original definition of P 5y 
obtaining the definition of the most general biharmonic fifth degree polynomial 

P 5 (x, y) = A 0 (i 5 - 5 xy 4 ) + A x (x 4 y - y 5 / 5) 

+A 2 (x 3 y 2 - xy 4 ) + A 3 (x 2 y 3 - y s /5) . (5.13) 

This function will be biharmonic for any values of the four independent constants 
Ao, A\, A 2 , A 2 . We can express this by stating that the biharmonic polynomial P 5 has 
four degrees of freedom. 

In general, the polynomial Q is of degree 4 less than P because the biharmonic 
equation is of degree 4. It follows that there are always four fewer constraint equations 
than there are coefficients in the original polynomial P and hence that they can 
be satisfied leaving a polynomial with 4 degrees of freedom. However, the process 
degenerates if N<3. 

In view of the above discussion, it might seem appropriate to write an expression 
for the general polynomial of degree N in the form of equation (5.13) as a preliminary 
to the solution of polynomial problems in rectangular coordinates. However, as can 
be seen from equation (5.13), the resulting expressions are algebraically messy and 
this approach becomes unmanageable for problems of any complexity. Instead, it 
turns out to be more straightforward algebraically to define problems in terms of the 
simpler unconstrained polynomials like equation (5.2) and to impose the constraint 
equations at a later stage in the solution. 

5.1.1 Second and third degree polynomials 

We recall that the stress components are defined in terms of the stress function <j) 
through the relations 



d 2 <t> 

° xx ~ dtf 
d 2 <p 

° yy ~ dx 2 

d 2 <t> 

° xy dxdy 



(5.14) 

(5.15) 

(5.16) 



It follows that when the stress function is a polynomial of degree N in x,y, the 
stress components will be polynomials of degree (N - 2). In particular, constant and 
linear terms in <j> correspond to null stress fields (zero stress everywhere) and can be 
disregarded. 

The second degree polynomial 



<t> ~ A 0 x 2 + A\xy + A 2 y' z 



(5.17) 
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yields the stress components 

@xx — 2A 2 , (7 X y — A\ } @yy — 2j4o (5.18) 

and hence corresponds to the most general state of biaxial uniform stress. 

The third degree polynomial 

(j> = A 0 x 3 + A\X 2 y + A 2 xy 2 + A 3 y 3 (5.19) 

yields the stress components 

cf xx = 2 A 2 x + 6A 3 y ; o xy = -2A x x - 2 A 2 y ; a yy = 6zloz + 2A x y . (5.20) 

If we arbitrarily set A 0 ,A X ,A 2 = 0, the only remaining non-zero stress will be 

o xx = 6A 3 y , (5.21) 

which corresponds to a state of pure bending, when applied to the rectangular beam 
—a<x<a, -b<y<b, as shown in Figure 5.1. 




Figure 5.1: The rectangular beam in pure bending. 

The other terms in equation (5.19) correspond to a more general state of bending. 
For example, the constant 2lo describes bending of the beam by tractions o yy applied 
to the boundaries y = ±b, whilst the terms involving shear stresses a xy could be 
obtained by describing a general state of biaxial bending with reference to a Cartesian 
coordinate system which is not aligned with the axes of the beam. 

The above solutions are of course very elementary, but we should remember that, 
in contrast to the Mechanics of Materials solutions for simple bending, they are 
obtained without making any simplifying assumptions about the stress fields. For ex- 
ample, we have not assumed that plane sections remain plane, nor have we demanded 
that the beam be long in comparison with its depth. Thus, the present section could 
be taken as verifying the exactness of the Mechanics of Materials solutions for uniform 
stress and simple bending, as applied to a rectangular beam. 




5.2. 
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5.2 Rectangular beam problems 

5.2.1 Bending of a beam by an end load 

Figure 5.2 shows a rectangular beam, 0 <x<a,-b<y<b, subjected to a transverse 
force, F at the end x = 0, and built-in at the end x = a, the horizontal boundaries 
•y = ±b being traction free. The boundary conditions for this problem are most 
naturally written in the form 



®xy 


= 0 ; 


•O 

-H 

II 


(5.22) 


a yy 


= 0 ; 


-o 

-H 

II 


(5.23) 


o xx 


= 0 ; 


x = 0 


(5.24) 


/ rT x yd y 
J —b 


= F 


; x = 0 . 


(5.25) 




Figure 5.2: Cantilever with an end load. 

The boundary condition (5.25) is imposed in the weak form, which means that 
the value of the traction is not specified at each point on the boundary — only the 
force resultant is specified. In general, we shall find that problems for the rectangular 
beam have finite polynomial solutions when the boundary conditions on the ends are 
stated in the weak form, but that the strong (i.e. pointwise) boundary condition can 
only be satisfied on all the boundaries by an infinite series or transform solution. This 
issue is further discussed in Chapter 6. 

Mechanics of Materials considerations suggest that the bending moment in this 
problem will vary with x and hence that the stress component a xx will have a leading 
term proportional to xy. This in turn suggests a fourth degree polynomial term xy 3 in 
the stress function (f>. Our procedure is therefore to start with the trial stress function 

<t> = <W , (5.26) 
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examine the corresponding tractions on the boundaries and then seek a corrective so- 
lution which, when superposed on equation (5.26), yields the solution to the problem. 
Substituting (5.26) into (5.14-5.16), we obtain 



&xx — QC\xy 


(5.27) 


& xy * == 3^x2/ 


(5.28) 


II 

0 


(5.29) 



from which we note that the boundary conditions (5.23, 5.24) are satisfied identically, 
but that (5.22) is not satisfied, since (5.28) implies the existence of an unwanted 
uniform shear traction - ZC\b 2 on both of the edges y = ±b. This unwanted traction 
can be removed by superposing an appropriate uniform shear stress, through the 
additional stress function term C 2 xy. Thus, if we define 

= Cixy 3 + C 2 xy , (5.30) 



equations (5.27, 5.29) remain unchanged, whilst (5.28) is modified to 

o xy = -3Ciy 2 - C 2 . (5.31) 



The boundary condition (5.22) can now be satisfied if we choose C 2 to satisfy the 
equation 

C 2 = -3Ci6 2 , (5.32) 

so that 

a xy = 3C 1 (5 2 - y 2 ) . (5.33) 

The constant C\ can be determined by substituting (5.33) into the remaining 
boundary condition (5.25), with the result 




The final stress field is therefore defined through the stress function 

F(xy 3 - 3 b 2 xy) 

V ' 

the corresponding stress components being 

3 Fxy 
2 b 3 

3F(b 2 - y 2 ) 

46 3 

0 . 



® XX — 

a xy = 

a yy = 



(5.34) 

(5.35) 

(5.36) 

(5.37) 

(5.38) 



The solution of this problem is given in the Mathematica and Maple files ‘S521’. 
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We note that no boundary conditions have been specified on the built-in end, 
x = a. In the weak form, these would be 



J &xx dy — 0 ; x = a 


(5.39) 


f a xy dy = F ; x = a 
J —b 


(5.40) 


J Oxxydy = Fa ; x = a . 


(5.41) 



However, if conditions (5.22-5.25) are satisfied, (5.39-5.41) are merely equivalent 
to the condition that the whole beam be in equilibrium. Now the Airy stress function 
is so defined that whatever stress function is used, the corresponding stress field will 
satisfy equilibrium in the local sense of equations (2.5). Furthermore, if every particle 
of a body is separately in equilibrium, it follows that the whole body will also be in 
equilibrium. It is therefore not necessary to enforce equations (5.39-5.41), since if we 
were to check them, we should necessarily find that they are satisfied identically. 

5.2.2 Higher order polynomials — a general strategy 

In the previous section, we developed the solution by trial and error, starting from 
the leading term whose form was dictated by equilibrium considerations. A more 
general technique is to identify the highest order polynomial term from equilibrium 
considerations and then write down the most general polynomial of that degree and 
below. The constant multipliers on the various terms are then obtained by imposing 
boundary conditions and biharmonic constraint equations. 

The only objection to this procedure is that it involves a lot of algebra. For 
example, in the problem of §5.2.1, we would have to write down the most general 
polynomial of degree 4 and below, which involves 12 separate terms even when we 
exclude the linear and constant terms as being null. However, this is not a serious 
difficulty if we are using Maple or Mathematica, so we shall first develop the steps 
needed for this general strategy. Shortcuts which would reduce the complexity of the 
algebra in a manual calculation will then be discussed in §5.2.3. 

Order of the polynomial 

Suppose we have a normal traction on the surface y = b varying with x n . In Mechanics 
of Materials terms, this corresponds to a distributed load proportional to x" and 
elementary equilibrium considerations show that the bending moment can then be 
expected to contain a term proportional to x n+2 . This in turn implies a bending stress 
a xx proportional to x n+2 y and a term in the stress function proportional to x n+2 y 3 — 
i.e. a term of polynomial order (n+5). A corresponding argument for shear tractions 
proportional to x m shows that we require a polynomial order of at least (m+4). 

A sufficient polynomial order can therefore be selected by the following procedure:- 
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(i) Identify the highest order polynomial term n in the normal tractions a yy on the 
surfaces y = ±b. 

(ii) Identify the highest order polynomial term m in the shear tractions a yx on the 
surfaces y = ±6. 

(iii) Use a polynomialfor 0 including all polynomial terms oforder max(m+4, n + 5) 
and below, but excluding constant and linear terms. 

In the special case where both surfaces y — ±b are traction-free, it is sufficient to use 
a polynomial of 4th degree and below (as in §5.2.1). 

Solution procedure 

Once an appropriate polynomial has been identified for 0, we proceed as follows:- 

(i) Substitute 0 into the biharmonic equation (5.1), leading to a set of constraint 
equations, as in §5.1. 

(ii) Substitute 0 into equations (5.14-5.16), to obtain the stress components as 
functions of x, y. 

(iii) Substitute the equations defining the boundaries (e.g. x = 0, y = b, y = -b in the 
problem of §5.2.1) into appropriate 1 stress components, to obtain the tractions 
on each boundary. 

(iv) For the longer boundaries (where strong boundary conditions will be enforced), 
sort the resulting expressions into powers of a: or y and equate coefficients with 
the corresponding expression for the prescribed tractions. 

(v) For the shorter boundaries, substitute the tractions into the appropriate weak 
boundary conditions, obtaining three further independent algebraic equations. 

The equations so obtained will generally not all be linearly independent, but they 
will be sufficient to determine all the coefficients uniquely. The solvers in Maple and 
Mathematic a can handle this redundancy. 

Example 

We illustrate this procedure with the example of Figure 5.3, in which a rectangular 
beam -a<x<a, — b<y<b is loaded by a uniform compressive normal traction p on 
y — b and simply supported at the ends. 

'Recall from §1.1.1 that the only stress components that act on (e.g.) y = 6 are those which 
contain y as one of the suffices. 
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Figure 5.3: Simply supported beam with a uniform load. 

The boundary conditions on the surfaces y — ±b can be written 

(5.42) 

(5.43) 

(5.44) 



Oy X 


~ 0 ; 


y = ±b 


a yy 


= ~P 


\ y = b 


°yy 


= 0 ; 


II 

1 



These boundary conditions are to be satisfied in the strong sense. To complete the 
problem definition, we shall require three linearly independent weak boundary con- 
ditions on one or both of the ends i = ±a. We might use symmetry and equilibrium 
to argue that the load will be equally shared between the supports, leading to the 
conditions 2 



F x {a) = _ 


rb 

/ a xx {a,y)dy = 0 

' —b 


(5.45) 




rb 

/ Oxy{a,y)dy = pa 

' —b 


(5.46) 


M(a) = / 


b 

°xx {a,y)ydy - 0 


(5.47) 



on the end x = a. As explained in §5.2.1, we do not need to enforce the additional 
three weak conditions on x — —a. 

The normal traction is uniform — i.e. it varies with a: 0 (n = 0), so the above 
criterion demands a polynomial of order (n + 5) = 5. We therefore write 

(f> = CiX 2 + C^xy +Czy 2 + C4X 3 + CsX 2 y+Cexy 2 + C7?/ 3 + Cgi 4 + C s x 3 y + C\oX 2 y 2 
■+C , nXi/ 3 +Ci2j/ 4 + Ci3X b +CuX i y+C\$x i y 2 +C\GX i y' i +Cnxy* +Ci&y 5 . (5.48) 

‘It is not necessary to use symmetry arguments to obtain three linearly independent weak con- 
ditions. Since the beam is simply supported, we know that 



Af(tt) = J a xx (a,y)ydy = 0 ; M(~a) 




-a,y)ydy = 0 



F x (a) = / o xx (a,y)dy = 0 . 
J-b 



It is easy to verify that these conditions lead to the same solution as (5.45-5.47). 
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This is a long expression, but remember we only have to type it in once to the 
computer file. We can cut and paste the expression in the solution of subsequent 
problems (and the reader can indeed cut and paste from the web file ‘polynomial’). 
Substituting (5.48) into the biharmonic equation (5.1), we obtain 

(120(7,3 + 24(7,5 + 24Cn)x + (24(7,4 + 24(7, 8 + 120(7, 8 )2/ + (24(7 8 + 8(7, „ + 24(7,2) = 0 

(5.49) 

and this must be zero for all x, y leading to the three constraint equations 

120Ci 3 + 24Cis + 24C ir = 0 (5.50) 

24Cu + 24Ci6 + 120C 18 = 0 (5.51) 

24(78 + 8(7,0 + 24(7,2 = 0. (5.52) 



The stresses are obtained by substituting (5.48) into (5.14-5.16) with the result 
o xx = 2C3 + 2C 3 x + 6C77/ + 2C\qx 2 + 6(7,, xy + 12(7, 2 J/ 2 + 2(7, 5a: 3 + &C\§x 2 y 



+l2C u xy 2 + 20 C l8 y 3 (5.53) 

a xy = -Ci - 2C 5 x - 2C 6 y - 3C 9 x 2 - 4C, 0 a;j/ - 3C,,j/ 2 - 4 C^x 3 - 6 C i5 x 2 y 

-6Ci 6 xy 2 - 4 Cny 3 (5.54) 

Oyy = 2(7, + 6 C 4 X + 2C*,y + 12Cgx 2 + 6C 9 xy + 2Cioy 2 + 20(7,3X^ + 12(7, \x 2 y 
+6C l5 xy 2 + 2C l6 y 3 . (5.55) 

The tractions on y = b are therefore 

Oy X — — 4(7,42: 3 — (3(7 9 + 6C, 5 6)a; 2 — (2C5 + 4(7,o6 + &Cigb 2 )x 

—(C2 + 2C56 + 3(7, ,6 2 + 4(7, 76 s ) (5.56) 

Oyy = 20C n x 3 + {12C S + \2C u b)x 2 + {6C 4 + 6Cgb + 6C l5 b 2 )x 

+(2(7, + 2C56 + +2(7, 0 5 2 + 2(7, 6 6 3 ) (5.57) 

and these must satisfy equations (5.42, 5.43) for all x, giving 

4(7 H = 0 (5.58) 

3(7 9 + 6(7, 5 6 = 0 (5.59) 

2 C 5 + 4(7,o6 + 6 C 156 2 = 0 (5.60) 

C 2 + 2Cgb + 3C,,6 2 + ACnb 3 = 0 (5.61) 

20(7,3 = 0 (5.62) 

12(7 8 + 12(7,46 = 0 (5.63) 

6 C 4 + 6(7 9 6 + 6 C, s 6 2 = 0 (5.64) 

2(7, + 2C 5 b + 2(7, 0 6 2 + 2(7, 6 6 3 = -p . (5.65) 
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A similar procedure for the edge y — ~b yields the additional equations 



3C 9 -6C l5 b = 0 


(5.66) 


2 C 5 — 4Cio£> 4- 60igb 2 = 0 


(5.67) 


C 2 - 2C 6 b + 3C n b 2 - 4C l7 b 3 = 0 


(5.68) 


12 C g - 12 C 14 6 = 0 


(5.69) 


6 C 4 — 6Cgb 4- 6 C 156 2 = 0 


(5.70) 


2 Ci - 2C s b + 2Cio & 2 - 2C 16 6 3 = 0 . 


(5.71) 



On x = a, we have 

a xx ~ 2 C 3 4- 2Cea + 6C7J/ + 2Ciofl 2 4- 6 Cnay + l2C\2y 3 4- 2 Ci 5 a 3 + 6Ci 9 a 2 y 

+l2C„ay 2 4- 20Ci8j/ 3 (5.72) 

a xy = — C 2 — 2C$a — 2 C§y — 3 C 9 a 2 — 4Ci 9 ay — 3 Cny 2 — 4Ci 4 a 3 — &C\^a 2 y 

-6 C l6 ay 2 - 4 C 17 y 3 . (5.73) 

Substituting into the weak conditions (5.45-5.47) and evaluating the integrals, we 
obtain the three additional equations 

4C36 + 4C$ab 4- 4CioQ, 2 b 4- 8C126 3 4- 4C' 15 a 3 6 4- 8C\ 7 ab 3 = 0 (5.74) 

-2C 2 b - 4C 5 ab - 6C 9 a 2 b - 2Cnb 3 - 8C u a 3 b - 4C l6 ab 3 = pa (5.75) 

4C 7 b 3 4- 4C u ab 3 + 4C l6 a 2 b 3 4- 8C lg 6 5 = 0. (5.76) 

Finally, we solve equations (5.50-5.52, 5.58-5.71, 5.74-5.76) for the unknown con- 
stants Ci, ..., Ci8 and substitute back into (5.48), obtaining 

<t> = — ^(5x 2 t/ 3 - y 3 - 15 b 2 x 2 y - 5 a 2 y 3 4- 2 b 2 y 3 - 1 0b 3 x 2 ) . (5.77) 

4Utr 

The corresponding stress field is 

**» = ^(I5x 2 y-10y 3 -15a 2 y + 6b 2 y) (5.78) 

° xy = (5J9) 

= £;(y 3 -3b 2 y-2b 3 ). (5.80) 

The reader is encouraged to run the Maple or Mathematica files ‘S522’, which 
contain the above solution procedure. Notice that most of the algebraic operations 
are generated by quite simple and repetitive commands. These will be essentially 
similar for any polynomial problem involving rectangular coordinates, so it is a simple 
matter to modify the program to cover other cases. 
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5.2.3 Manual solutions — symmetry considerations 

If the solution is to be obtained manually, the complexity of the algebra makes the 
process time consuming and increases the likelihood of errors. Fortunately, the com- 
plexity can be reduced by utilizing the natural symmetry of the rectangular beam. In 
many problems, the loading has some symmetry which can be exploited in limiting 
the number of independent polynomial terms and even when this is not the case, 
some saving of complexity can be achieved by representing the loading as the sum of 
symmetric and antisymmetric parts. We shall illustrate this procedure by repeating 
the solution of the problem of Figure 5.3. 

The problem is symmetrical about the mid-point of the beam and hence, taking 
the origin there, we deduce that the resulting stress function will contain only even 
powers of x. This immediately reduces the number of terms in the general stress 
function to 10. 

The beam is also symmetrical about the axis y — 0, but the loading is not. We 
therefore decompose the problem into the two sub-problems illustrated in Figure 
5.4(o,6). 



Oyy = -p/2 



Gyy - -p/2 



pa 




Oyy -p/2 



Oyy “ ~p/2 



(a) (b) 

Figure 5.4: Decomposition of the problem into (a) antisymmetric and ( b ) symmetric 
parts. 

The problem in Figure 5.4(a) is antisymmetric in y and hence requires a stress 
function with only odd powers of y , whereas that of Figure 5.4(b) is symmetric and 
requires only even powers. In fact, the problem of Figure 5.4(b) clearly has the trivial 
solution corresponding to uniform uniaxial compression, o yy = —p/2, the appropriate 
stress function being cj)=—px 2 / 4. 

For the problem of Figure 5.4(a), the most general fifth degree polynomial which 
is even in x and odd in y can be written 

<t> = C$x 2 y + Cyy 3 + CuX A y + Ci&x 2 y 3 4- Cigj/ 5 , 



(5.81) 
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which has just five degrees of freedom. We have used the same notation for the re- 
maining constants as in (5.48) to aid in comparing the two solutions. The appropriate 
boundary conditions for this sub-problem are 



O xy — 0 , 


y = ±b 


(5.82) 


«ci 

II 

-H 

to i*« 


; y = ±b 


(5.83) 


J b Vxxdy = 0 ; 


x = ±a 


(5.84) 


f b 

J b o xx ydy = 0 ; 


x = ±a . 


(5.85) 



Notice that, in view of the symmetry, it is only necessary to satisfy these conditions 
on one of each pair of edges (e.g. on y = b, x = a). For the same reason, we do not have 
to impose a condition on the vertical force at x = ±a, since the symmetry demands 
that the forces be equal at the two ends and the total force must be 2 pa to preserve 
global equilibrium, this being guaranteed by the use of the Airy stress function, as in 
the problem of §5.2.1. 

It is usually better strategy to start a manual solution with the strong bound- 
ary conditions (equations (5.82, 5.83)), and in particular with those conditions that 
are homogeneous (in this case equation (5.82)), since these will often require that 
one or more of the constants be zero, reducing the complexity of subsequent steps. 
Substituting (5.81) into (5.15, 5.16), we find 

o xy = -2C$x - 4C 14 :r 3 — QCiexy 2 (5.86) 

o yy = 2C 5 y + 12C u x 2 y + 2C l6 y 3 ■ (5.87) 

Thus, condition (5.82) requires that 

4C u x 3 + (2 C 5 + 6Ci 6 6 2 )a; = 0 ; for all x (5.88) 

and this condition is satisfied if and only if 

Ci 4 = 0 and 2C S + 6C 16 5 2 = 0 . (5.89) 

A similar procedure with equation (5.87) and boundary condition (5.83) gives the 
additional equation 

2C 5 b + 2Ci 6 b 3 = . (5.90) 

Equations (5.89, 5.90) have the solution 

We next determine C\& from the condition that the function (j> is biharmonic, 
obtaining 

(24C 14 + 24C i6 + 120C lg )y = 0 (5.92) 
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and hence 




(5.93) 



from (5.89, 5.91, 5.92). 

It remains to satisfy the two weak boundary conditions (5.84, 5.85) on the ends 
x = ±a. The first of these is satisfied identically in view of the antisymmetry of the 
stress field and the second gives the equation 



4 C 7 b 3 + 4C 16 g 2 5 3 + 8 C 18 /> 5 = 0 , 



(5.94) 



which, with equations (5.91, 5.93), serves to determine the remaining constant, 



C 7 



p(2b 2 — 5a 2 ) 
406 5 



(5.95) 



The final solution of the complete problem (the sum of that for Figures 5.4(a) and 
(b)) is therefore obtained from the stress function 

<j) = (5z 2 y 3 - y 5 - 15 b 2 x 2 y - 5 a 2 y 3 + 2 b 7 y 3 - 10fc 3 x 2 ) , (5.96) 

as in the ‘computer solution’ (5.77), and the stresses are therefore given by (5.78-5.80) 
as before. 



5.3 Fourier series and transform solutions 

Polynomial solutions can, in principle, be extended to more general loading of the 
beam edges, as long as the tractions are capable of a power series expansion. However, 
the practical use of this method is limited by the algebraic complexity encountered for 
higher order polynomials and by the fact that many important traction distributions 
do not have convergent power series representations. 

A more useful method in such cases is to build up a general solution by components 
of Fourier form. For example, if we write 

<t> = f{y) cos(Xx) or <b = f(y) sin(Ax) , (5.97) 

substitution in the biharmonic equation (5.1) shows that f(y) must have the general 
form 

f(y) = (A + By)e Xy + (C + Dy)e~** , (5.98) 

where A, B, C, D are arbitrary constants. Alternatively, by defining new arbitrary 
constants A',B',C',D' through the relations A — A! + C', B — B' + D' , C = A' - 
C, D = B'-D', we can group the exponentials into hyperbolic functions, obtaining 
the equivalent form 



f{y) = (A' + B'y) cosh (Ay) + {C + D'y ) sinh(Ay) . 



(5.99) 
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The hyperbolic form enables us to take advantage of any symmetry about y = 0, 
since cosh(Ay), y sinh(Ay) are even functions of y and sinh(Ay), ycosh(Ay) are odd 
functions. 

More general biharmonic stress functions can be constructed by superposition of 
terms like (5.98, 5.99), leading to Fourier series expansions for the tractions on the 
surfaces y = ±b. The theory of Fourier series can then be used to determine the 
coefficients in the series, using strong boundary conditions on y = ±6. Quite general 
traction distributions can be expanded in this way, so Fourier series solutions provide 
a methodology applicable to any problem for the rectangular bar. 

5.3.1 Choice of form 

The stresses due to the stress function 4> = f(y) cos(Ax) are 

o xx = -A 2 f{y) cos(Ax) ; o xy = A f'{y) sin(Ax) ; a yy = f"[y) cos(Ax) (5.100) 

and the tractions on the edge x = a are 

o xx {a, y) = -A 2 /(y) cos(Aa) ; o xy (a , y) = A/'(y) sin(Aa) . (5.101) 

It follows that we can satisfy homogeneous boundary conditions on one (but not both) 
of these tractions in the strong sense, by restricting the Fourier series to specific values 
of A. In equation (5.101), the choice A = 7rn/a will give a xy = 0 on x = ±a, whilst 
A = (2 n — l)7r/2 a will give cr xx =0 on x=±a, where n is any integer. 

Example 

We illustrate this technique by considering the rectangular beam — a < x < a, -b < 
y<b , simply supported at x — ±a and loaded by compressive normal tractions Pi(x) 
on the upper edge y = b and P 2 (x) on y = — b — i.e. 



Oxy = o ; y = ±b 


(5.102) 


Oyy = ~Pl( x ) ; y = b 


(5.103) 


= -pi{x) ; y = -6 


(5.104) 


o xx = 0 ; x = ±a . 


(5.105) 



Notice that we have replaced the weak conditions (5.84, 5.85) by the strong condition 
(5.105). As in §5.2.2, it is not necessary to enforce the remaining weak conditions 
(those involving the vertical forces on x = ±o), since these will be identically satisfied 
by virtue of the equilibrium condition. 

The algebraic complexity of the problem will be reduced if we use the geometric 
symmetry of the beam to decompose the problem into four sub-problems. For this 
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purpose, we define 



M x ) = fi(~ x ) = 


^{Pl{x) + Pi(-x) + p 2 (x) + p 2 (-x)} 


(5.106) 


f2{x) = -/ 2 (“X) = 


- Pl{~ x ) + P2(x) -Pii-x)} 


(5.107) 


M x ) = h(~ x ) = 


j{pi(x) +pi(-x) -p 2 (x) p 2 ( x) } 


(5.108) 


M x ) = -fii-x) = 


^{p,(x) -Pt(-x) -P 2 (x) +p 2 (-x)} 


(5.109) 

(5.110) 



and hence 

Pi(z) = /i(z) + / 2 (x)+/ 3 (x)+/ 4 (x) ; P 2 {x) = Moc) + f 2 {x)-f 3 (x)-U{x) . (5.111) 
The boundary conditions now take the form 



o xy = 0 ; y = ±b (5.112) 

a n = ~M X ) - M x ) - M x ) - M x ) ; y = * (5.i 13) 

= -fi(x) - / 2 (x) + f 3 (x) + / 4 (x) ; y = -b (5.114) 

a xx = 0 ; x = ±a (5.115) 



and each of the functions /i, / 2 , / 3 , / 4 defines a separate problem with either symmetry 
or antisymmetry about the x- and y- axes. We shall here restrict attention to the 
loading defined by the function / 3 (a;), which is symmetric in x and antisymmetric in 
y. The boundary conditions of this sub-problem are 



a xy = 0 ; y = ±b (5.116) 

°yy = =F/3(z) i V=±b (5.117) 

a xx = 0 ; x = ±a . (5.118) 

The problem of equations (5.116-5.118) is even in x and odd in y, so we use a 
cosine series in x with only the odd terms from the hyperbolic form (5.99) — i.e. 

OO 

<t>= X)M„ycosh(A n y) + £ n sinh(A„t/)}cos(A n i) , (5.119) 

n=l 

where A n ,B n are arbitrary constants. The strong condition (5.118) on x — ±a can 
then be satisfied in every term by choosing 



x (2 n — l)7r 
" = 2a 



(5.120) 
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The corresponding stresses are 

OO 

o xx = 5Z{24 n A n sinh(A n t/)+A n A^2/cosh(A n y) + 5 n A^sinh(A n t/)}cos(A n x) (5.121) 

n= 1 
oo 

o xy = y{A n X n cosh(X n y)+A n X 2 n ysinh(X n y)+B n X 2 n cosh{X n y)}sm{X n x) (5.122) 

n=l 

oo 

°yy ~ 53{ j4 nA^ycosh(A„y) + B n A* sinh(A„y)} cos(A n x) (5.123) 

n=l 

and hence the boundary conditions (5.116, 5.117) on y = ±b require that 

OO 

EM n A n cosh(X n b) + A n Xnbsmh(X n b)+B n Xl cosh(A n 5)} sin(A n j) = 0 (5.124) 

n=l 

oo 

EMnA^6 cosh(A„6) + Z? n A^sinh(A n 6)} cos(A„x) = / 3 (x) .(5.125) 

n=l 



To invert the series, we multiply (5.125) by cos(A m x) and integrate from —a to a, 
obtaining 



00 * a 

T. / {^l n A^6cosh(A n 6) 4- B n A^sinh(A„6)}cos(A n x)cos(A m x)dx 

= [ f 3 (x) cos(X m x)dx . (5.126) 

J —a 



The integrals on the left hand side are all zero except for the case m = n and hence, 
evaluating the integrals, we find 

{A m \ 2 m bcosh{X m b) + B m X 2 m s\nh{X m b)}a= f f 3 (x) cos{X m x)dx . (5.127) 

J —a 

The homogeneous equation (5.124) is clearly satisfied if 

A m X m cosh(A m 6) 4- >l m A^6sinh(A m i) + B m X 2 m cosh(A m 6) = 0 . (5.128) 



Solving (5.127, 5.128) for we have 

cosh(A m 5) 



B m = 



J f 3 (x) cos(X m x)dx 

m i\ u • , ,, nT / /s(*) cos(X m x)dx 

X 2 n o.{X m b - sinh(A m 6)cosh(A m 6)} J- a 



A m a{A m 6 - sinh(A m 6) cosh(A m 6)} 
(cosh(A m 6) + A m 6sinh(A m 6)) 



(5.129) 



where X„, is given by (5.120). The stresses are then recovered by substitution into 
equations (5.121-5.123). 

The corresponding solutions for the functions /i ,/ 2 ,/4 are obtained in a sim- 
ilar way, but using a sine series for the odd functions / 2,/4 and the even terms 
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ysinh(Ay), cosh(Ay) in <p for /^/j. The complete solution is then obtained by super- 
posing the solutions of the four sub-problems. 

The Fourier series method is particularly useful in problems where the traction 
distribution on the long edges has no power series expansion, typically because of 
discontinuities in the loading. For example, suppose the beam is loaded only by a 
concentrated compressive force F on the upper edge at x = 0, corresponding to the 
loading p\(x) = F5(x), Pi(x) = 0. For the symmetric/antisymmetric sub-problem 
considered above, we then have 



/*(*) = ^ (5.130) 

from (5.108) and the integral in equations (5.129) is therefore 

J f 3 (x) cos{X m x)dx = j . (5.131) 

This solution satisfies the end condition on a xx in the strong sense, but the con- 
dition on a xy only in the weak sense. In other words, the tractions a xy on the ends 
add up to the forces required to maintain equilibrium, but we have no control over 
the exact distribution of these tractions. This represents an improvement over the 
polynomial solution of §5.2.3, where weak conditions were used for both end tractions, 
so we might be tempted to use a Fourier series even for problems with continuous 
polynomial loading. However, this improvement is made at the cost of an infinite 
series solution. If the series were truncated at a finite value of n, errors would be 
obtained particularly near the ends or any discontinuities in the loading. 

5.3.2 Fourier transforms 

If the beam is infinite or semi-infinite (a— >oo), the series (5.119) must be replaced 
by the integral representation 

roo 

<t>{x,y)= / f{\,y)cos{\x)d\ , (5.132) 

Jo 

where 

f{X,y) - A(A)ycosh(Ay) + B(A) sinh(Ay) . (5.133) 

Equation (5.132) is introduced here as a generalization of (5.97) by superposition, but 
(f)(x,y) is in fact the Fourier cosine transform of /(A, y), the corresponding inversion 
being 

2 r°° 

f(X,y) = - <f>(x, y) cos(Ax)cte . (5.134) 

7C Jo 

The boundary conditions on y = ±b will also lead to Fourier integrals, which can 
be inverted in the same way to determine the functions A(X)>B( A). For a definitive 




PROBLEMS 



67 



treatment of the Fourier transform method, the reader is referred to the treatise by 
Sneddon 3 . Extensive tables of Fourier transforms and their inversions are given by 
Erdelyi 4 . The cosine transform (5.132) will lead to a symmetric solution. For more 
general loading, the complex exponential transform can be used. 

It is worth remarking on the way in which the series and transform solutions are 
natural generalizations of the elementary solution (5.97). One of the most powerful 
techniques in Elasticity — and indeed in any physical theory characterized by linear 
partial differential equations — is to seek a simple form of solution (often in separated- 
variable form) containing a parameter which can take a range of values. A more 
general solution can then be developed by superposing arbitrary multiples of the 
solution with different values of the parameter. 

For example, if a particular solution can be written symbolically as <j>= 
where A is a parameter, we can develop a general series form 

OO 

4>{x,y) = Y, A 'f( x >y> x i) (5.135) 

i - 0 

or an integral form 

4>{x,y)= f A{\)f(x,y,\)d\ . (5.136) 

Ja 

The series form will naturally arise if there is a discrete set of eigenvalues, A j for 
which f(x, y, A;) satisfies some of the boundary conditions of the problem. Additional 
examples of this kind will be found in §§6.2, 11.2. In this case, the series (5.135) is most 
properly seen as an eigenfunction expansion. Integral forms arise most commonly 
(but not exclusively) in problems involving infinite or semi-infinite domains (see, for 
example, §§11.3, 27.2.2.). 

Any particular solution containing a parameter can be used in this way and, since 
transforms are commonly named after their originators, the reader desirous of instant 
immortality might like to explore some of those which have not so far been used. Of 
course, the usefulness of the resulting solution depends upon its completeness — i.e. 
its capacity to represent all stress fields of a given class — and upon the ease with 
which the transform can be inverted. 

PROBLEMS 

1. The beam —b<y<b, 0 < x < L, is built-in at the end x = 0 and loaded by a 
uniform shear traction o xy = S on the upper edge, y = b, the remaining edges, x = 
L, y = — 6 being traction-free. Find a suitable stress function and the corresponding 
stress components for this problem, using the weak boundary conditions on x = L. 

3 I.N.Sneddon, Fourier Transforms, McGraw-Hill, New York, 1951. 

4 A.Erdelyi, ed., Tables of Integral Transforms, Bateman Manuscript Project, California Institute 
of Technology, Vol.l, McGraw-Hill, New York, 1954. 
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2. The beam —b <y < b, -L < x < L is simply supported at the ends x = ±L and 
loaded by a shear traction o xy = Sx/L on the lower edge, y = -b , the upper edge being 
traction-free. Find a suitable stress function and the corresponding stress components 
for this problem, using the weak boundary conditions on x = ±L. 

3. The beam —b < y <b, 0 < x < L, is built-in at the end x = L and loaded by a 
linearly-varying compressive normal traction p(x) = Sx/L on the upper edge, y = b, 
the remaining edges, x = 0, y = — b being traction-free. Find a suitable stress function 
and the corresponding stress components for this problem, using the weak boundary 
conditions on x = 0. 

4. The beam -b < y < b, - L < x < L is simply supported at the ends x = ±L and 
loaded by a compressive normal traction 

V(X) = 5 cos (g) 

on the upper edge, y =■ b, the lower edge being traction-free. Find a suitable stress 
function and the corresponding stress components for this problem. 

5. The beam —b<y<b, 0 < x < L, is built-in at the end x = L and loaded by a 
compressive normal traction 

PW = Ssin (g) 

on the upper edge, y — b , the remaining edges, x = 0, y = —b being traction-free. Use 
a combination of the stress function (5.97) and an appropriate polynomial to find the 
stress components for this problem, using the weak boundary conditions on x = 0. 

6. A large plate defined by y>0 is subjected to a sinusoidally varying load 

Oyy = S sin Ax ; a xy = 0 



at its plane edge y = 0. 

Find the complete stress field in the plate and hence estimate the depth y at which 
the amplitude of the variation in o yy has fallen to 10% of S. 

Hint: You might find it easier initially to consider the case of the layer 0 <y<h, 
with y = h traction-free, and then let /;. — > oo. 

7. The beam — a < x < a, -b < y < b is loaded by a uniform compressive traction p in 
the central region —a/2 < x < a/2 of both of the edges y — ±b, as shown in Figure 
5.5. The remaining edges are traction-free. Use a Fourier series with the appropriate 
symmetries to obtain a solution for the stress field, using the weak condition on a xy 
on the edges x=±a and the strong form of all the remaining boundary conditions. 
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8. Use a Fourier series to solve the problem of Figure 5.4(a) in §5.2.3. Choose the 
terms in the series so as to satisfy the condition a xx (±a,y) = 0 in the strong sense. 

If you are solving this problem in Maple or Mathematica, compare the solution 
with that of §5.2.3 by making a contour plot of the difference between the truncated 
Fourier series stress function and the polynomial stress function 

<t> = -~p (5xV -y 5 - 15 b 2 x 2 y - 5a 2 y 3 + 2 b 2 y 3 ) . 

Examine the effect of taking different numbers of terms in the series. 

9. The large plate y > 0 is loaded at its remote boundaries so as to produce a state 
of uniform tensile stress 

Oxx — ^ i ® xy — ty — d i 

the boundary y = 0 being traction-free. We now wish to determine the perturbation 
in this simple state of stress that will be produced if the traction-free boundary had 
a slight waviness, defined by the line 

y = ecos(A:r) , 

where Ac <C 1. To solve this problem 
(i) Start with the stress function 



C ,,2 

<fi = ~y + f(y) cos(Ax) 

and determine f(y) if the function is to be biharmonic. 

(ii) The perturbation will be localized near y = 0, so select only those terms in /( y) 
that decay as y— Kx>. 
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(iii) Find the stress components and use the stress transformation equations to de- 
termine the tractions on the wavy boundary. Notice that the inclination of 
the wavy surface to the plane y = 0 will be everywhere small if Xe 1 and 
hence the trigonometric functions involving this angle can be approximated 
using sin(x)«:r, cos(i)wl,ar<l. 

(iv) Choose the free constants in f(y) to satisfy the traction-free boundary condition 
on the wavy surface. 

(v) Determine the maximum tensile stress and hence the stress concentration factor 
as a function of Xe. 




Chapter 6 

END EFFECTS 



The solution of §5.2.2 must be deemed approximate insofar as the boundary conditions 
on the ends x = ±a of the rectangular beam are satisfied only in the weak sense of 
force resultants, through equation (5.45-5.47). In general, if a rectangular beam is 
loaded by tractions of finite polynomial form, a finite polynomial solution can be 
obtained which satisfies the boundary conditions in the strong (i.e. pointwise) sense 
on two edges and in the weak sense on the other two edges. 

The error involved in such an approximation corresponds to the solution of a 
corrective problem in which the beam is loaded by the difference between the actual 
tractions applied and those implied by the approximation. These tractions will of 
course be confined to the edges on which the weak boundary conditions were applied 
and will be self-equilibrated, since the weak conditions imply that the tractions in 
the approximate solution have the same force resultants as the actual tractions. 

For the particular problem of §5.2.2, we note that the stress field of equations 
(5.78-5.80) satisfies the boundary conditions on the edges y = ±b, but that there is a 
self-equilibrated normal traction 

axx ~W& tttfy-by 3 ) ( 6 . 1 ) 

on the ends x = ±o, which must be removed by superposing a corrective solution if 
we wish to satisfy the boundary conditions of Figure 5.3 in the strong sense. 



6.1 Decaying solutions 

In view of Saint- Venant’s theorem, we anticipate that the stresses in the corrective 
solution will decay as we move away from the edges where the self-equilibrated trac- 
tions are applied. The decay rate is likely to be related to the width of the loaded 
region and hence we anticipate that the stresses in the corrective solution will be 
significant only in two regions near the ends, of linear dimensions comparable to the 
width of the beam. These regions are, shown shaded in Figure 6.1. It follows that 
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the solution of §5.2.2 will be a good approximation in the unshaded region in Figure 

6.1. 



imuiiii 



pa 




\pa 



Figure 6.1: Regions of the beam influenced by end effects. 

It also follows that the corrective solutions for the two ends are uncoupled, since 
the corrective field for the end x = —a has decayed to negligible proportions before 
we reach the end x = +a and vice versa. This implies that, so far as the left end 
is concerned, the corrective solution is essentially identical to that which would be 
required in the semi-infinite beam, x > —a. For the rest of this discussion, we shall 
therefore simplify the statement of the problem by considering the corrective solution 
for the left end only and shifting the origin to —a, so that the semi-infinite beam 
under consideration is defined by i>0, —b<y<b. 

It is also now clear why we chose to satisfy the strong boundary conditions on the 
long edges, y = ±b. If instead we had imposed strong conditions on x = ±a and weak 
conditions on y = ±6, the shaded regions would have overlapped and there would be 
no region in which the finite polynomial solution was a good approximation to the 
stresses in the beam. It also follows that the approximation will only be useful when 
the beam has an aspect ratio significantly different from unity — i.e. b/a^l. 



6.2 The corrective solution 

We recall that the stress function, <f>, for the corrective solution must (i) satisfy the 
biharmonic equation (5.1), (ii) have zero tractions on the boundaries y = ±b and (iii) 
have prescribed non-zero tractions (such as those defined by equation (6.1)) on the 
end(s). We cannot generally expect to find a solution to satisfy all these conditions in 
closed form and hence we seek a series or transform (integral) solution as suggested 
in §5.3. However, the solution will be simpler if we can find a class of solutions, all of 
which satisfy some of the conditions. We can then write down a more general solution 
as a superposition of such solutions and choose the coefficients so as to satisfy the 
remaining condition(s). 

Since the traction boundary condition on the end will vary from problem to prob- 
lem, it is convenient to seek solutions which satisfy (i) and (ii) — i.e. biharmonic 
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functions which define zero tractions on the surfaces y = ±b. The traction-free condi- 
tion requires 

&yy — &yx — 0 , y — (6.2) 



and hence 



d 2 <j> 
dx 2 



d 2 <t> 

dxdy 



= 0 ; y = ±b. 



(6.3) 



6.2.1 Separated-variable solutions 

One way to obtain functions satisfying conditions (6.3) is to write them in the 
separated-variable form 

4> = f{x)g(y) , (6.4) 

in which case, (6.3) will be satisfied for all x, provided that 

g(y) = g'iy) = 0 ; y = ±b. (6.5) 



Notice that the final corrective solution cannot be expected to be of separated- 
variable form, but it is possible that it can be represented as the sum of such terms. 
If the functions (6.4) are to be biharmonic, we must have 



9 dx* dx 2 dy 2 J dy* ' 



( 6 . 6 ) 



and this equation must be satisfied for all values of a;, y. Now, if we consider the subset 
of points x.c. where c is a constant, it is clear that f(x) must satisfy an equation of 
the form 

- 4 ^ + B 0 +c/ = 0 ' (67) 
where A, B, C, are constants, and hence f(x) must consist of exponential terms such 
as f(x) = exp(Ai). Similar considerations apply to the function g(y). Notice inci- 
dentally that X might be complex or imaginary, giving sinusoidal functions, and there 
are also degenerate cases where C and/or B = 0 in which case /( x) could also be a 
polynomial of degree 3 or below. 

Since we are seeking to represent a field which decays with x, we select terms of 
the form 

d> = 9(y)e~ Xx , (6.8) 

in which case, (6.6) reduces to 



+ 2 vf| + v s = 0 . 



(6.9) 



d 4 g 

-r— r ~r 4* , 
dy* dy 2 

which is a fourth order ordinary differential equation for g(y) with general solution 



g{y) = (4i + A 2 y) cos(Xy) + (A 3 + A t y) sin {Xy) . (6.10) 
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6.2.2 The eigenvalue problem 

The arbitrary constants A[, A 2 , A3, A 4 are determined from the boundary conditions 
(6.2), which in view of (6.5) lead to the four simultaneous equations 

(A\ + A^b) cos(A 6 ) + (A 3 + A 4 y) sin(A 6 ) = 0 (6.11) 

(ylj - A 2 b) cos(Afc) - (A 3 - A 4 y) sin(A 6 ) = 0 (6.12) 

( A 2 + A3X H- A 4 Xb) cos(A b) — (y4i A + A 2 A 6 — A 4 ) sin(A 6 ) = 0 (6.13) 

(A? + A3X — A 4 Xb) cos(A5) 4- (AiA - A 2 \b - At)(sin A b) — 0 . (6.14) 

This set of equations is homogeneous and will generally have only the trivial solu- 
tion A\ = A 2 — ^3 = A) =0. However, there are some eigenvalues for the exponential 
decay rate, A, for which the determinant of coefficients is singular and the solution is 
non-trivial. 

A more convenient form of the equations can be obtained by taking sums and dif- 
ferences in pairs — i.e. by constructing the equations (6.11+6.12), (6.11-6.12), (6.13+ 
6.14), (6.13 - 6.14), which after rearrangement and cancellation of non-zero factors 
yields the set 





A\ cos(A b) + A|6sin(A6) 


= 0 


(6.15) 


A] Xsin(Xb) 


— yl 4 {sin(A6) -+ A6(cos(A6)} 


= 0 


(6.16) 




A 2 bcos(Xb) + A3 sin(A6) 


= 0 


(6.17) 


A 2 {cos(A6) 


— A6sin(A6)} + A3\c.os(\b) 


= 0 . 


(6.18) 



What we have done here is to use the symmetry of the system to partition the 
matrix of coefficients. The terms A|Cos(Ay), Aiysm(Xy) are symmetric, whereas 
A 2 ycos(\y), A3sin(Ay) are antisymmetric. The boundary conditions are also sym- 
metric and hence the symmetric and antisymmetric terms must separately satisfy 
them. 

We conclude that the set of equations (6.15-6.18) has two sets of eigenvalues, for 
one of which the resulting eigenfunction is symmetric and the other antisymmetric. 
The symmetric eigenvalues are obtained by eliminating A U A 4 from (6.15, 6.16) with 
the result 

sin(2A6) + 2A6 = 0 , (6.19) 

whilst the antisymmetric eigenvalues are obtained in the same way from (6.17, 6.18) 
with the result 

sin(2A6) - 2A6 = 0 . (6.20) 

Figure 6.2 demonstrates graphically that the only real solution of equations (6.19, 
6.20) is the trivial case A = 0 (which in fact corresponds to the non-decaying solutions 
in which a force or moment resultant is applied at the end and transmitted along the 
beam). 
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Figure 6.2: Graphical solution of equations (6.19, 6.20). 

However, there is a denumerably infinite set of non-trivial complex solutions, cor- 
responding to stress fields which oscillate whilst decaying along the beam. These 
solutions are fairly easy to find by writing Xb — c + td, separating real and imaginary 
parts in the complex equation, and solving the resulting two simultaneous equations 
for the real numbers, c, d, using a suitable numerical algorithm. 

Once the eigenvalues and the related eigenfunctions have been determined, the 
next stage is to establish a more general solution of the form 

OO 

i=i 

where gi(y) is the eigenfunction corresponding to the eigenvalue and C{ defines a 
set of as yet undetermined coefficients. 

The final step is to choose the constants C, so as to satisfy the prescribed boundary 
condition on the end x = 0. This raises a question of completeness. It is fairly clear 
that (i) if we truncate the series (6.21), we can develop an approximate solution to 
the problem and (ii) that the accuracy of the solution (however defined) can always 
be improved by taking more terms. In particular cases, this is easily established 
numerically. However, it is more challenging to prove that the eigenfunction expansion 
is complete in the sense that any prescribed self-equilibrated traction on x = 0 can 
be described to within an arbitrarily small accuracy by taking a sufficient number of 
terms in the series, though experience with other eigenfunction expansions (e.g. with 
expansion of elastodynamic states of a structure in terms of normal modes) suggests 
that this will always be true. In fact, although the analysis described in this section 
has been known since the investigations by Papkovich and Fadle in the early 1940s, 
the formal proof of completeness was only completed by Gregory 1 in 1980. 

*R.D. Gregory, Green's functions, bi-linear forms and completeness of the eigenfunctions for the 
elastostatic strip and wedge, J. Elasticity. Vol. 9 (1979), 283-309; R.D. Gregory, The semi-infinite 
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It is worth noting that, as in many related problems, the eigenfunctions oscillate 
in y with increasing frequency as A* increases and in fact every time we increase i by 
1, an extra zero appears in the function <ji(y) in the range 0<y<b. Thus, there is a 
certain similarity to the process of approximating functions by Fourier series and in 
particular, the residual error in case of truncation will always cross zero once more 
than the last eigenfunction included. 

This is also helpful in that it enables us to estimate the decay rate of the first 
excluded term. We see from equation (6.10) that the distance between zeros in y in 
any of the separate terms would be {it /\ r), where \r is the real part of A. Itfollows 
that over a corresponding distance in the x direction, the field would decay by the 
factor exp(— 7t) = 0.0432. This suggests that we might estimate the decay rate of 
the end field by noting the distance between zeros in the corresponding tractions 2 . 
For example, the traction of equation (6.1) has zeros at y = 0,3b/5, corresponding to 
A*=(5ir/36) = 5.23/6. 

An alternative way of estimating the decay rate is to note that the decay rate for 
the various terms in (6.21) increases with i and hence as x increases, the leading term 
will tend to predominate. The tractions of (6.1) are antisymmetric and hence the 
leading term corresponds to the real part of the first eigenvalue of equation (6.20), 
which is found numerically to be A/j6 = 3.7. 

Either way, we can conclude that the error associated with the end tractions in 
the approximate solution of §5.2.2 has decayed to around e -4 , i.e. to about 2% of the 
values at the end, within a distance b of the end. Thus the region affected by the end 
condition — the shaded region in Figure 6.1 — is quite small. 



(a) 




(b) 
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Figure 6.3: Leading term in self-equilibrated tractions for (a) symmetric loading 
and (b) antisymmetric loading. 

strip x > 0, — 1 < y < 1; completeness of the Papkovich-Fadle eigenfunctions when0 IX (O, J/), 0yy (0, y) 
are prescribed, J. Elasticity, Vol. 10 (1980), 57-80; R.D. Gregory, The traction boundary value prob- 
lem for the elastostatic semi-infinite strip; existence of solution and completeness of the Papkovich- 
Fadle eigenfunctions, J. Elasticity, Vol. 10 (1980), 295-327. These papers also include extensive 
references to earlier investigations of the problem. 

2 This assumes that the wavelength of the tractions is the same as that of <j>, which of course is 
an approximation, since neither function is purely sinusoidal. 
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For problems which are symmetric in y, the leading self-equilibrated term is likely 
to have the form of Figure 6.3(a), which has a longer wavelength than the correspond- 
ing antisymmetric form, 6.3(b). The end effects in symmetric problems therefore decay 
more slowly and this is confirmed by the fact that the real part of the first eigenvalue 
of the symmetric equation (6.19) is only A«6 = 2.1. 



6.3 Other Saint- Venant problems 

The general strategy used in this chapter can be used in other curvilinear coordinate 
systems to correct the errors incurred by imposing the weak boundary conditions on 
appropriate edges. The essential steps are:- 

(i) Define a coordinate system (£, rj) such that the boundaries on which the strong 
conditions are applied are of the form, ? 7 =constant. 

(ii) Find a class of separated-variable biharmonic functions containing a parameter 
(A in the above case). 

(iii) Set up a system of four homogeneous equations for the coefficients of each 
function, based on the four traction-free boundary conditions for the corrective 
solution on the edges y =constant. 

(iv) Find the eigenvalues of the parameter for which the system has a non-trivial 
solution and the corresponding eigenfunctions, which are then used as the terms 
in an eigenfunction expansion to define a general form for the corrective field. 

(v) Determine the coefficients in the eigenfunction expansion from the prescribed 
inhomogeneous boundary conditions on the end £ constant. 

PROBLEMS 

1. Show that if £ = x+iy and sin(C)-C = 0. where x,y are real variables, then 

f(x) = cosx\/x 2 - sin 2 ! + sin x ln(sin x) — sin x In + \]x 2 - sin 2 x j = 0 . 

Using Maple or Mathematica to plot the function f(x), find the first six roots of this 
equation and hence determine the first six values of Xnb for the antisymmetric mode. 

2. Devise a method similar to that outlined in Problem 6. 1 to determine the first six 
values of A nb for the symmetric mode. 
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3. A displacement function representation for plane strain problems can be developed 3 
in terms of two harmonic functions <£, u in the form 

„ d(p du> d<t> dw , n i . 

2 "“- = ai + ; 2 "“« = % + ^ ( w 



@xx 



d 2 (j) d 2 u> du 

dx 2 + y dx 2 dy 



Oxy — 



d 2 (j) 

dxdy 



d 2 u> du> 

+ V dxdy ( 2 ^dx 



a yy 



d 2 (f> d 2 u> . 



v) 



dui 
dy ' 



Use this representation to formulate the eigenvalue problem of the long strip x > 
0 ,—b<y<b whose edges y = ±b are both bonded to a rigid body. Find the eigenvalue 
equation for symmetric and antisymmetric modes and comment on the expected decay 
rates for loading of the strip on the end x = 0. 



4. Use the displacement function representation of Problem 6.3 to formulate the 
eigenvalue problem for the long strip x > 0, -b < y < b whose edges y = ±6 are 
in frictionless contact with a rigid body (so that the normal displacement is zero, 
but the frictional (tangential) traction is zero). Find the eigenvalue equation for 
symmetric and antisymmetric modes and comment on the expected decay rates for 
loading of the strip on the end x = 0. 



5. Use the displacement function representation of Problem 6.3 to formulate the 
eigenvalue problem for the long strip x > 0, -b < y < b which is bonded to a rigid 
surface at y = — b , the other long edge y — b being traction-free. Notice that this 
problem is not symmetrical, so the problem will not partition into symmetric and 
antisymmetric modes. 



3 This solution is a two-dimensional version of a three-dimensional solution developed in Chapter 
19 and tabulated in Table 19.1 as solutions A and B. 
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A body force is defined as one which acts directly on the interior particles of the body, 
rather than on the boundary. Since the interior of the body is not accessible, it follows 
necessarily that body forces can only be produced by some kind of physical process 
which acts ‘at a distance’. The commonest examples are forces due to gravity and 
magnetic or electrostatic attraction. In addition, we can formulate elasticity problems 
for accelerating bodies in terms of body forces, using D’Alembert’s principle (see 
§7.2.2 below). 



7.1 Stress function formulation 



We noted in §4.4 that the Airy stress function formulation satisfies the equilibrium 
equations if and only if the body forces are identically zero, but the method can be 
extended to the case of non-zero body forces provided the latter can be expressed as 
the gradient of a scalar potential, V. 

We adopt the new definitions 



a xx 

a yy 



Try — 





d 2 <f> 

dxdy 



(7.1) 

(7.2) 

(7.3) 



in which case the two-dimensional equilibrium equations will be satisfied provided 



dV dV 

Vx ~ dx ' Py ~ dy ’ 



(7.4) 



i.e. 



p= -VK. 



(7.5) 
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Notice that the body force potential V appears in the definitions of the two normal 
stress components a xx ,a yy (7.1, 7.2), but not in the shear stress a xy (7.3). Thus, the 
modification in these equations is equivalent to the addition of a biaxial hydrostatic 
tension of magnitude V, which is of course invariant under coordinate transformation 
(see §1.1.3). 



7.1.1 Conservative vector fields 



If a force field is capable of being represented as the gradient of a scalar potential, 
as in equation (7.5), it is referred to as conservative. This terminology arises from 
gravitational theory, since if we move a particle around in a gravitational field and if 
the force on the particle varies with position, the principle of conservation of energy 
demands that the work done in moving the particle from point A to point B should be 
path-independent — or equivalently, the work done in moving it around a closed path 
should be zero. If this were not the case, we could choose a direction for the particle 
to move around the path which would release energy and hence have an inexhaustible 
source of energy. 

If all such integrals are path-independent, we can use the work done in bringing 
a particle from infinity to a given point as the definition of a unique local potential. 
Then, by equating the work done in an infinitesimal motion to the corresponding 
change in potential energy, we can show that the local force is proportional to the 
gradient of the potential, thus demonstrating that a conservative force field must 
be capable of a representation like (7.5). Conversely, if a given force field can be 
represented in this form, we can show by integration that the work done in moving a 
particle from A to B is proportional to V(A)-V{B) and is therefore path-independent. 

Not all body force fields are conservative and hence the formulation of §7.1 is not 
sufficiently general for all problems. However, we shall show below that most of the 
important problems involving body forces can be so treated. 

We can develop a condition for a vector field to be conservative in the same way 
as we developed the compatibility conditions for strains. We argue that the two 
independent body force components p x ,p y are defined in terms of a single scalar 
potential V and hence we can obtain a constraint equation on p x ,p y by eliminating 
V between equations (7.4) with the result 



dp y dp x 
dx dy 



(7.6) 



In three dimensions there are three equations like (7.6), from which we conclude 
that a vector field pis conservative if and only if curl p = 0. Another name for such 
fields is irrotational, since we note that if we replace p by u, the conditions like (7.6) 
are equivalent to the statement that the rotation u> is identically zero (cf equation 
(1.40)). 

If the body force field satisfies equation (7.6), the corresponding potential can be 
recovered by partial integration. We shall illustrate this procedure in §7.2 below. 
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7.1.2 The compatibility condition 



We demonstrated in §4.4.3 that, in the absence of body forces, the compatibility 
condition reduces to the requirement that the Airy stress function <f> be biharmonic. 
This condition is modified when body forces are present. 

We follow the same procedure as in §4.4.3, but use equations (7. 1-7.3) in place 
of (4.1). Substituting into the compatibility equation in terms of stresses (4.4), we 
obtain 



o>> d 2 V d*<j> d 2 V dU 

dy 4 dy 2 dx 2 dy 2 dy 2 dx 2 dy 2 

d*<t> d 2 V d 4 (t> _ SV 

+ dx* dx 2 dx 2 dy 2 dx 2 



(7.7) 



i.e. 



VV = -( l-t/)V 2 K. 



(7.8) 



Methods of obtaining suitable functions which satisfy this equation will be dis- 
cussed in §7.3 below. 



7.2 Particular cases 

It is worth noting that the vast majority of mechanical engineering components are 
loaded principally by boundary tractions rather than body forces. Of course, most 
components are subject to gravity loading, but the boundary loads are generally 
so much larger that gravity can be neglected. This is less true for civil engineering 
structures such as buildings, where the self-weight of the structure may be much larger 
than the weight of the contents or wind loads, but even in this case it is important 
to distinguish between the gravity loading on the individual component and that 
transmitted to the component by way of boundary tractions. 

It might be instructive at this point for the reader to draw free-body diagrams for a 
few common engineering components and identify the sources and relative magnitudes 
of the forces acting upon them. There are really comparatively few ways of applying 
a load to a body. By far the commonest is to push against it with another body — 
in other words to apply the load by contact. This is why contact problems occupy 
a central place in elasticity theory 1 . Significant loads may also be applied by fluid 
pressure as in the case of turbine blades or aircraft wings. Notice that it is fairly easy 
to apply a compressive normal traction to a boundary, but much harder to apply 
tension or shear. 

This preamble might be taken as a justification for not studying the subject of 
body forces at all, but there are a few applications in which they are of critical 
importance, most notably those dynamic problems in which large accelerations occur. 



'See Chapters 12, 26. 
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We shall develop expressions for the body force potential for some important cases in 
the following sections. 

7.2.1 Gravitational loading 

The simplest type of body force loading is that due to a gravitational field. If the 
problem is to remain two-dimensional, the direction of the gravitational force must 
lie in the xy-plane and we can choose it to be in the negative y-direction without loss 
of generality. The magnitude of the force will be pg per unit volume, where p is the 
density of the material, so in the notation of §2.1 we have 

Px — 0 ; Py = ~pg • ( 7 - 9 ) 

This force field clearly satisfies condition (7.6) and is therefore conservative and 
by inspection we note that it can be derived from the body force potential 

v = pgy ■ (7.10) 

It also follows that V 2 V = 0 and hence the stress function, <j> for problems involving 
gravitational loading is biharmonic, from (7.6). 

7.2.2 Inertia forces 

It might be argued that D’Alembert achieved immortality simply by moving a term 
from one side of an equation to the other, since D'Alembert’s principle consists merely 
of writing Newton’s second law of motion in the form F — mo=0 and treating the 
term —ma as a fictitious force in order to reduce the dynamic problem to one in 
statics. 

This simple process enables us to formulate elasticity problems for accelerating 
bodies as elastostatic body force problems, the corresponding body forces being 

Px = ~pa x ; Py~~pay, (7.11) 

where a x ,a y are the local components of acceleration, which may of course vary with 
position through the body. 

7.2.3 Quasi-static problems 

If the body were rigid, the accelerations of equation (7.11) would be restricted to those 
associated with rigid-body translation and rotation, but in a deformable body, the 
distance between two points can change, giving rise to additional, stress-dependent 
terms in the accelerations. 

These two effects give qualitatively distinct behaviour, both mathematically and 
physically. In the former case, the accelerations will generally be defined a priori 
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from the kinematics of the problem, which therefore reduces to an elasticity problem 
with body forces. We shall refer to such problems as quasi-static. 

By contrast, when the accelerations associated with deformations are important, 
they are not known a priori, since the stresses producing the deformations are them- 
selves part of the solution. In this case the kinematic and elastic problems are cou- 
pled and must be solved together. The resulting equations are those governing the 
propagation of elastic waves through a solid body and their study is known as Elas- 
todynamics. 

In this chapter, we shall restrict attention to quasi-static problems. As a practical 
point, we note that the characteristic time scale of elastodynamic problems is very 
short. For example it generally takes only a very short time for an elastic wave to 
traverse a solid. If the applied loads are applied gradually in comparison with this 
time scale, the quasi-static assumption generally gives good results. A measure of the 
success of this approximation is that it works quite well even for the case of elastic 
impact between bodies, which may have a duration of the order of a few milliseconds 2 . 

7.2.4 Rigid-body kinematics 

The most general acceleration for a rigid body in the plane involves arbitrary trans- 
lation and rotation. We choose a coordinate system fixed in the body and suppose 
that, at some instant, the origin has velocity Vo and acceleration a 0 and the body is 
rotating in the clockwise sense with absolute angular velocity Q and angular acceler- 
ation Q. The instantaneous acceleration of the point (r, 6) relative to the origin can 
then be written 

a T = -fl 2 r ; ag = -tlr . ( 712 ) 

Transforming these results into the x, ^-coordinate system and adding the accel- 
eration of the origin, we obtain the components of acceleration of the point ( x,y ) 
as 

a x = oq x - Cl 2 x + tly ( 7 . 13 ) 

a y = a 0y — Cl 2 y — tlx ( 7 . 14 ) 

and hence the corresponding body force field is 

Vx = -p(a 0x -tfx + tly) ( 7 . 15 ) 

p y = -p(a,oy - Sl 2 y - tlx) . ( 7 . 16 ) 

The astute reader will notice that the case of gravitational loading can be recovered 
as a special case of these results by writing ao y = g and setting all the other terms to 

2 For more information about elastodynamic problems, the reader is referred to the classical 
texts of J.D.Achenbach, Wave Propagation in Elastic Solids, North Holland, Amsterdam, (1973) 
and A.C.Eringen and E.S.fjuhubi, Elastodynamics, Academic Press, New York, (1975). For a more 
detailed discussion of the impact of elastic bodies, see K.L.Johnson, Contact Mechanics, Cambridge 
University Press, Cambridge, (1985), §11.4. 




84 



CHAPTER 7. BODY FORCES 



zero. In fact, a reasonable interpretation of the gravitational force is as a D’Alembert 
force consequent on resisting the gravitational acceleration. Notice that if a body is 
in free fall — i.e. if it is accelerating freely in a gravitational field and not rotating — 
there is no body force and hence no internal stress unless the boundaries are loaded. 

Substitution of (7.15, 7.16) into (7.6) shows that the inertia forces due to rigid- 
body accelerations are conservative if and only if 0 = 0 — i.e. if the angular velocity is 
constant. We shall determine the body force potential for this special case. Methods 
of treating the problem with non-zero angular acceleration are discussed in §7.4 below. 
From equations (7.4, 7.15, 7.16) with 0 = 0 we have 

p{a 0x - 0 2 x) (7.17) 

p(a 0y - 0 2 y) , (7.18) 

and hence, on partial integration of (7.17) 

V = p - ^0 2 £ 2 j + h{y ) , (7.19) 

where h(y) is an arbitrary function of i/only. Substituting this result into (7.18) we 
obtain the ordinary differential equation 

^ = (7-20) 

for h(y), which has the general solution 



dV 

dx 

dV 

dy 



h{y) = p (a 0y y - 4- C , (7.21) 

where C is an arbitrary constant which can be taken to be zero without loss of 
generality, since we are only seeking a particular potential function V. 

The final expression for V is therefore 

V = P (aoxX + a 0y y - ^0 2 (i 2 + y 2 )) . (7.22) 

The reader might like to try this procedure on a set of body forces which do not 
satisfy the condition (7.6). It will be found that the right hand side of the ordinary 
differential equation like (7.20) then contains terms which depend on x and hence 
this equation cannot be solved for h(y). 
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7.3 Solution for the stress function 

Once the body force potential V has been determined, the next step is to find a 
suitable function </>, which satisfies the compatibility condition (7.8) and which de- 
fines stresses through equations (7. 1-7.3) satisfying the boundary conditions of the 
problem. There are broadly speaking two ways of doing this. One is to choose some 
suitable form (such as a polynomial) without regard to equation (7.8) and then satisfy 
the constraint conditions resulting from (7.8) in the same step as those arising from 
the boundary conditions. The other is to seek a general solution of the inhomoge- 
neous equation (7.8) and then determine the resulting arbitrary constants from the 
boundary conditions. 

7.3.1 The rotating rectangular beam 

As an illustration of the first method, we consider the problem of the rectangular 
beam -a < x < a, —b<y<b, rotating about the origin at constant angular velocity 
f2, all the boundaries being traction-free (see Figure 7.1). 




Figure 7.1: The rotating rectangular bar. 

The body force potential for this problem is obtained from equation (7.22) as 

V = -±fi&(x 2 + t?), (7-23) 

and hence the stress function must satisfy the equation 

V 4 <p = 2p(l — i/)Q 2 , (7.24) 

from (7.8, 7.23). 

The geometry suggests a formulation in Cartesian coordinates and equation (7.24) 
leads us to expect a polynomial of degree 4 in x,y. We also note that V is even in 
both x and y and that the boundary conditions are homogeneous, so we propose the 
candidate stress function 



<j) = A x x 1 + A 2 x 2 y 2 + A 3 y 4 + A 4 X 1 + A h y 2 



(7.25) 




CHAPTER 7. BODY FORCES 



which contains all the terms of degree 4 and below with the required symmetry. 

The constants will be determined from equation (7.24) and from the 

boundary conditions 



<JyX 


= 0 ; 


y = ±6 


(7.26) 


°yy 


= 0 ; 


y = ±b 


(7.27) 


1 ^ o xx dy 


= 0 ; 


x = ±a , 


(7.28) 



where we have applied the weak boundary conditions only on x = ±a. Note also that 
the other two weak boundary conditions — that there should be no moment and no 
shear force on the ends — are satisfied identically in view of the symmetry of the 
problem about t/ = 0. 

As in the problem of §5.2.3, it is algebraically simpler to start by satisfying the 
strong boundary conditions (7.26, 7.27). Substituting (7.25) into (7. 1-7.3), we obtain 

a xx = 2A 2 x 2 + 12A 3 y 2 + 2A b - ^pil 2 {x 2 + y 2 ) (7.29) 

o yy = 12 A lX 2 + 2 A 2 y 2 + 2 A i -^pQ 2 (x 2 + y 2 ) (7.30) 

Oy X = -4 A 2 xy . (7-31) 

It follows that conditions (7.26, 7.27) will be satisfied for all x if and only if 

Ai = pi} 2 / 24 ; A 2 = 0 ; A 4 = pQ 2 b 2 / 4 . (7.32) 

The constant A 3 can now be determined by substituting (7.25) into (7.24), with 
the result 

24(Ai + A 3 ) = 2pft 2 (l - v ) , (7.33) 

which is the inhomogeneous equivalent of the constraint equations (see §5.1) and 
hence 

A 3 = pft 2 (l - 2i/)/24 , (7.34) 

using (7.32). 

Finally, we determine the remaining constant A 5 by substituting (7.29) into the 
weak boundary condition (7.28) and evaluating the integral, with the result 

* = '® 2 (t + t)- (7 ' 35) 




(7.36) 

(7.37) 

(7.38) 
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from equations (7.29-7.31). 

Notice that the boundary conditions on the ends x = ±a agree with those of the 
physical problem except for the second term in a xx , which represents a symmetric 
self-equilibrated traction. From §6.2.2, we anticipate that the error due to this dis- 
agreement will be confined to regions near the ends of length comparable with the 
half-length, b. 

7.3.2 Solution of the governing equation 

Equation (7.24) is an inhomogeneous partial differential equation — i.e. it has a 
known function of x, y on the right hand side — and it can be solved in the same way 
as an inhomogeneous ordinary differential equation, by first finding a particular solu- 
tion of the equation and then superposing the general solution of the corresponding 
homogeneous equation. 

In this context, a particular solution is any function <p that satisfies (7.24). It 
contains no arbitrary constants. The generality in the general solution comes from 
arbitrary constants in the homogeneous solution. Furthermore, the homogeneous 
solution is the solution of equation (7.24) modified to make the right hand side zero 
— i.e. 

V> = 0 (7.39) 

Thus, the homogeneous solution is a general biharmonic function and we can 
summarize the solution method as containing the three steps:- 

(i) finding any function <f> which satisfies (7.24); 

(ii) superposing a sufficiently general biharmonic function (which therefore contains 
several arbitrary constants); 

(iii) choosing the arbitrary constants to satisfy the boundary conditions. 

In the problem of the preceding section, the particular solution would be any 
fourth degree polynomial satisfying (7.24) and the homogeneous solution, the most 
general fourth degree biharmonic function with the appropriate symmetry. 

Notice that the particular solution (i) is itself a solution of a different physical 
problem — one in which the body forces are correctly represented, but the correct 
boundary conditions are not (usually) satisfied. Thus, the function of the homoge- 
neous solution is to introduce additional degrees of freedom to enable us to satisfy 
the boundary conditions. 

Physical superposition 

It is often helpful to think of this superposition process in a physical rather than a 
mathematical sense. In other words, we devise a related problem in which the body 
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forces are the same as in the real problem, but where the boundary conditions are 
simpler. For example, if a beam is subjected to gravitational loading, a simple physical 
‘particular’ solution would correspond to the problem in which the beam is resting on 
a rigid foundation and hence the stress field is one-dimensional. To complete the real 
problem, we would then have to superpose the solution of a corrective problem with 
tractions equal to the difference between the required tractions and those implied by 
the particular solution, but with no body force, since this has already been taken into 
account in the particular solution. 

One advantage of this way of thinking is that it is not restricted to problems in 
which the body force can be represented by a potential. We can therefore use it to 
solve the problem of the rotationally accelerating beam in the next section. 



7.4 Rotational acceleration 

We saw in §7.2.4 that the body force potential cannot be used in problems where 
the angular acceleration is non-zero. In this section, we shall generate a particular 
solution for this problem and then generalize it, using the results without body force 
from Chapter 5. 

7.4.1 The circular disk 

Consider the rotationally symmetric problem of Figure 7.2. A solid circular disk, 
radius a, is initially at rest (fi = 0) and at time t = 0 it is caused to accelerate in 
a clockwise direction with angular acceleration by tractions uniformly distributed 
around the edger =a. Note that we could determine the magnitude of these tractions 
by writing the equation of motion for the disk, but it will not be necessary to do this 
— the result will emerge from the analysis of the stress field. 



“ O>0 




Figure 7.2: The disk with rotational acceleration. 

At any given time, the body forces and hence the resulting stress field will be the 
sum of two parts, one due to the instantaneous angular velocity and the other to the 
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angular acceleration. The former can be obtained using the body force potential and 
we therefore concentrate here on the contribution due to the angular acceleration. 
This is the only body force at the beginning of the process, so the following solution 
can also be regarded as the solution for the instant, ( = 0. 

The problem is clearly axisymmetric, so the stresses and displacements must de- 
pend on the radius r only, but it is also antisymmetric, since if the sign of the tractions 
were changed, the problem would become a mirror image of that illustrated. Now sup- 
pose some point in the disk had a non-zero outward radial displacement, u r . Changing 
the sign of the tractions would change the sign of this displacement — i.e. make it be 
directed inwards — but this is impossible if the new problem is to be a mirror image 
of the old. We can therefore conclude from symmetry that u r = 0 throughout the 
disk. In the same way, we can conclude, that the stress components o rr , ogg are zero 
everywhere. Incidentally, we also note that if the problem is conceived as being one 
of plane strain, symmetry demands that o zz be everywhere zero. It follows that this 
is one of those lucky problems in which plane stress and plane strain are the same 
and hence the plane stress assumption involves no approximation. 

We conclude that there is only one non-zero displacement component, ug, and 
one non-zero strain component, e r g. Thus, the number of strain and displacement 
components is equal and no non-trivial compatibility conditions can be obtained by 
eliminating the displacement components. (An alternative statement is that all the 
compatibility equations are satisfied identically, as can be verified by substitution — 
the compatibility equations in cylindrical polar coordinates are given by Saada 3 ). 

It follows that the only non- zero stress component, o T g can be determined from 
equilibrium considerations alone. Considering the equilibrium of a small element due 
to forces in the 0-direction and dropping terms which are zero due to the symmetry 
of the system, we obtain 

do r $ 2o r0 • 

— + — = -prtt , (7.40) 

which has the general solution 



o r e = - 



pi ir 2 
4 




(7.41) 



where A is an arbitrary constant which must be set to zero to retain continuity at 
the origin. Thus, the stress field in the disk of Figure 7.2 is 



pClr 2 

°r0 = — 

In particular, the traction at the surface r = a is 



-o r g{a) = 



ptla 2 



(7.42) 



(7.43) 



3 A.S. Saada, Elasticity, Pergamon Press, New York, (1973), §6.9. 
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and the applied moment about the axis of rotation is therefore 

M — — 27 ca 2 cr r o(a) = , (7.44) 

since the traction acts over a length 27ra and the moment arm is a. We know from 
elementary dynamics that M = I Cl, where I is the moment of inertia and we can 
therefore deduce from (7.44) that 



npa 4 
~ 2 ~ ’ 



(7.45) 



which is of course the correct expression for the moment of inertia of a solid disk of 
density p, radius a and unit thickness. Notice however that we were able to deduce 
the relation between M and without using the equations of rigid-body dynamics. 
Equation (7.40) ensures that every particle of the body obeys Newton’s law and this 
is sufficient to ensure that the complete body satisfies the equations of rigid-body 
dynamics. 



7.4.2 The rectangular bar 

We can use the stress field (7.42) as a particular solution for determining the stresses 
in a body of any shape due to angular acceleration. One way to think of this is to 
imagine cutting out the real body from an imaginary disk of sufficiently large radius. 
The stresses in the cut-out body will be the same as those in the disk, provided we 
arrange to apply tractions to the boundaries of the body that are equal to the stress 
components on that surface before the cut was made. These will not generally be 
the correct boundary tractions for the real problem, but we can correct the boundary 
tractions by superposing a homogeneous solution — i.e. a corrective solution for 
the actual body with prescribed boundary tractions (equal to the difference between 
those applied and those obtained in the disk solution) but no body forces (since these 
have already been taken care of in the disk solution.) 

As an illustration, we consider the rectangular bar, —a < x < a, —b < y < b, 
accelerated by two equal shear forces, F, applied at the ends x — ±a as shown in 
Figure 7.3, the other boundaries, y — ±b being traction free. 




Figure 7.3: The rectangular bar with rotational acceleration. 
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We first transform the disk solution into rectangular coordinates, using the trans- 
formation relations (1. 7-1.9), obtaining 



&XX 

°yy 

°xy 



- 2 cr r 0 sin 0 cost? = 
2 o r y sin 0 cos 0 = — 
(7 r fl(cos 2 0 - sin 2 0) 



ptlr 2 sin 0 cos 0 
2 

pClr 2 sin 0 cos 0 
2 



ptlxy 

2 

pClxy 

2 ~ 



pflr 2 (cos 2 0 - sin 2 0) 
2 



(7.46) 



(7.47) 



/3Q(x 2 - y 2 ) 
4 



. (7.48) 



This stress field is clearly odd in both x and y and involves normal tractions on 
y = ±b which vary linearly with x. The bending moment will therefore vary with 
a; 3 suggesting a stress function (f> with a leading term x 5 y — i.e. a sixth degree 
polynomial. 

The most general polynomial with the appropriate symmetry is 

<j) - Aix 5 y + A 2 x 3 y 3 + A 3 xy 5 + A t x 3 y + A b xy 3 + A e xy . (7.49) 



The stress components are the sum of those obtained from the homogeneous stress 
function (7.49) using the definitions (4.1) — remember the homogeneous solution here 
is one without body force — and those from the disk problem given by equations 
(7.46-7.48). We find 



°xx = l,P&xy + 6 A 2 x 3 y + 20A 3 xy 3 + 6A$xy 


(7.50) 


Oyy = — ^pClxy + 2QA\X 3 y + 6A 2 xy 3 + 6A t xy 


(7.51) 


o xy = ~^p(l(x 2 - y 2 ) - 5AiX 4 - 9A 2 x 2 y 2 


-5A 3 y 4 - ZAiX 2 - ZA b y 2 - A 6 . 


(7.52) 


The boundary conditions are 


Oy X = 0 ; y~±b 


(7.53) 


<*: 

II 

O 

II 

H- 

o- 


(7.54) 


J y<*xxdy = o ; x ~ ±a , 


(7.55) 



where we note that weak boundary conditions are imposed on the ends x = ±a. Since 
the solution is odd in y, only the moment and shear force conditions are non-trivial 
and the latter need not be explicitly imposed, since they will be satisfied by global 
equilibrium (as in the problem of §5.2.2). 

Conditions (7.53, 7.54) have to be satisfied for all x and hence the corresponding 
coefficients of all powers of x must be zero. It follows immediately that A\ = 0 and 
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the remaining conditions can be written 



— -pfi + 6A 2 b 2 + 6i4< 
~ptl - 9A 2 b 2 - 3A 4 
\pYlb 2 - 5A 3 b 4 - 3A 5 b 2 - A 6 



0 

0 

0 



(7.56) 

(7.57) 

(7.58) 



We get one additional condition from the requirement that <j> (equation (7.49)) be 
biharmonic 

72i4 a + 120^3 = 0 (7.59) 



(7.60) 



and another from the boundary condition (7.55), which with (7.50) yields 

— ptl + 2A 2 q 2 ~b 4 ^ 36 ^ -F 2A 3 — 0 . 

6 

The solution of these equations is routine, giving the stress function 

<p = —~^r(5x 2 y 3 -3xy 5 — l0b 2 x 3 y + llb 2 xy 3 -5a 2 xy 3 + l5a 2 b 2 xy-33b 4 xy) . (7.61) 
60cr 

The complete stress field, including the particular solution terms, is 

■ 2 3 (a 2 



@XX 



°yy ~ 






2 ~ x2 )A 

2 b 2 ) 



pQxy 



l-t 

b 2 



a xy = pti(b 2 -y 2 ) 



y 2 + a 2 — 3x 2 11 



4 b 2 



20 y ’ 



(7.62) 

(7.63) 

(7.64) 



from (7.50-7.52). 

Finally, we can determine the forces F on the ends by integrating the shear trac- 
tion, a xy over either end as 

F = - f a xy (a)dy = ~p£lb(a 2 + b 2 ) . (7.65) 

J —b O 

Maple and Mathematica solutions of this problem are given in the files ‘S742’. As 
in §7.4.1, the relation between the applied loading and the angular acceleration has 
been obtained without recourse to the equations of rigid-body dynamics. However, 
we note that the moment of inertia of the rectangular bar for rotation about the 
origin is I = 4pab(a 2 + b 2 )/3 and the applied moment is M = 2 Fa, so application of 
the equation M = /17 leads to the same expression for F as that found in (7.65). 




PROBLEMS 



93 



7.4.3 Weak boundary conditions and the equation of motion 

We saw in §5.2.1 that weak boundary conditions need only be applied at one end of 
a stationary rectangular bar, since the stress field defined in terms of the Airy stress 
function must involve tractions that maintain the body in equilibrium. Similarly, 
in dynamic problems involving prescribed accelerations, an appropriate set of weak 
boundary conditions can be omitted. For example, we did not have to specify the 
value of F in the problem of Figure 7.3. Instead, it’s value was calculated after the 
stress field had been determined and necessarily proved to be consistent with the 
equations of rigid-body dynamics. 

If a body is prevented from moving by a statically determinate support, it is 
natural to treat the support reactions as the ‘neglected’ weak boundary conditions. 
Thus, if the body is attached to a rigid support at one boundary, we apply no weak 
conditions at that boundary, as in §5.2.1. If the body is simply supported at a 
boundary, we impose the weak boundary condition that there be zero moment applied 
at the support, leading for example to the conditions stated in footnote 2 on Page 59. 

Similar considerations apply in a dynamic problem if the body is attached to 
a support which moves in such a way as to prescribe the acceleration of the body. 
However, we may also wish to solve problems in which specified non-equilibrated loads 
are applied to an unsupported body. For example, we may be asked to determine the 
stresses in the body of Figure 7.3 due to prescribed end loads F. One way to do this 
would be first to solve a rigid-body dynamics problem to determine the accelerations 
and then proceed as in §7.4.2. However, in view of the present discussion, a more 
natural approach would be to include the angular acceleration 0 as an unknown and 
use equation (7.65) to determine it in terms of F at the end of the solution procedure. 

More generally, if we have a two-dimensional body subjected to prescribed trac- 
tions on all edges, we could assume the most general accelerations (7.13, 7.14) and 
solve the body force problem, treating a nx , a 0y , Ll as if they were known. If strong 
boundary conditions are imposed on two opposite edges and weak boundary condi- 
tions on both the remaining edges, it will then be found that there are three extra 
conditions which serve to determine the unknown accelerations. 



PROBLEMS 

1. Every particle of an elastic body of density p experiences a force 

„ CSm 



directed towards the origin, where C is a constant, r is the distance from the origin 
and 8m is the mass of the particle. Find a body force potential V that satisfies these 
conditions. 
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2. Verify that the body force distribution 

Pi = Cy ; p y = -Cx 

is non-conservative, by substituting into equation (7.6). Use the technique of §7.2.4 
to attempt to construct a body force potential V for this case. Identify the step at 
which the procedure breaks down. 



3. To construct a particular solution for the stress components in plane strain due 
to a non-conservative body force distribution, it is proposed to start by representing 
the displacement components in the form 



1 dip 
Ux ~ 2/i dy 



1 dtp 
2p dx 



, Uj. — 0 . 



Use the strain-displacement equations (1.43) and Hooke's law (1.60) to find expres- 
sions for the stress components in terms of ip. Substitute these results into the equi- 
librium equations (2.2, 2.3) to find the governing equations for the stress function 
ip. 

What is the condition that must be satisfied by the body force distribution p if 
these equations are to have a solution? Show that this condition is satisfied if the 
body force distribution can be written in terms of a potential function W as 



_ dW 
dy 



Py = 



dW 

dx 



For the special case 



Ps — Cy ; p y = -Cx , 



find a particular solution for ip in the axisymmetric polynomial form 



ip = A(x 2 + y 2 ) n , 



where A is a constant and n is an appropriate integer power. Show that this solution 
can be used to obtain the stress components (7.46-7.48). Suggest ways in which 
this method might be adapted to give a particular solution for more general non- 
conservative body force distributions. 



4. If the elastic displacement u varies in time, there will generally be accelerations 
a = u and hence body forces p = —pu, from equation (7.11). Use this result and 
equation (2.17) to develop the general equation of linear Elastodynamics. 

Show that this equation is satisfied by a displacement field of the form 

u x = f(x - cit) ; u y = g(x - c 2 t) ; u t = 0 , 

where /, g are any functions and Ci,C 2 are two constants that depend on the material 
properties. Find the values ofci,C 2 and comment on the physical significance of this 
solution. 
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5. The beam —b<y<b, 0<x<L is built-in at the edge x = L and is loaded merely 
by its own weight, pg per unit volume. 

Find a solution for the stress field, using weak conditions on the end x = 0. 

6. One wall of a multistory building of height H is approximated as a thin plate 
—b<x<b, 0 < y < H. During an earthquake, ground motion causes the building to 
experience a uniform acceleration a in the x-direction. Find the resulting stresses in 
the wall if the material has density p and the edges x = ±6, y — H can be regarded as 
traction-free. 

7. A tall thin rectangular plate -a < x < a, —b<y<b(b^>a) is supported on 
the vertical edges x = ±a and loaded only by its own weight (density p). Find the 
stresses in the plate using weak boundary conditions on the horizontal edges y = ±b 
and assuming that the support tractions consist only of uniform shear. 

8. Figure 7.4(a) shows a triangular cantilever, defined by the boundaries y = 0, y — 
x tana and built-in at x = a. It is loaded by its own weight, pg per unit volume. Find 
a solution for the complete stress field and compare the maximum tensile bending 
stress with that predicted by the Mechanics of Materials theory. 

Would the maximum tensile stress be lower if the alternative configuration of 
Figure 1.4(b) were used? 





9. The thin rectangular plate —a<x<a,—b<y<b with a b rotates about the 
y-axis at constant angular velocity fi. All surfaces of the plate are traction-free. Find 
a solution for the stress field, using strong boundary conditions on the long edges 
y = ±b and weak boundary conditions on the ends x~±a. 

10. Solve Problem 7.9 for the case where a<&b. In this case you should use strong 
boundary conditions on x = ±a and weak boundary conditions on y = ±b. 
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11. A thin triangular plate bounded by the lines y = ±x tana, x = a rotates about the 
axis x = a at constant angular velocity fi. The inclined edges of the plate y = ±z tan a 
are traction-free. Find a solution for the stress field, using strong boundary conditions 
on the inclined edges (weak boundary conditions will then be implied on x = a). 

12. A thin triangular plate bounded by the lines y = ±x tana, x = a is initially at rest, 
but is subjected to tractions at the edge x — a causing an angular acceleration 0 about 
the perpendicular axis through the point (a, 0). The inclined surfaces y = ±x tana 
are traction-free. Find a solution for the instantaneous stress field, using the technique 
of §7.4.2. 




Chapter 8 

PROBLEMS IN POLAR 
COORDINATES 



Polar coordinates (r, 0 ) are particularly suited to problems in which the boundaries 
can be expressed in terms of equations like r = a,0 = a. This includes the stresses in 
a circular disk or around a circular hole, the curved bar with circular boundaries and 
the wedge, all of which will be discussed in this and subsequent chapters. 



8.1 Expressions for stress components 



We first have to transform the biharmonic equation (4.7) and the expressions for 
stress components (4.1) into polar coordinates, using the relations 



x = r cos 8 ; y = r sin 8 



( 8 . 1 ) 



r = 




6 = arctan 




( 8 . 2 ) 



The derivation is tedious, but routine. We first note by differentiation that 



d 


dr d 


dO d 


.3 


sin# d 


dx 


dx dr dx dO 


= COS9 dr ~ 


r dO 


d 


dr d 


dO d 


a d 


cos# d 


dy 


dy dr + dy dO 


— sin#— + 
dr 


r dd 



It follows that 

dx 2 



. d sin 9 d 

cm9 fr-—re 



n d 

COS 0- 

dr 



2 a d 2 . 2 f 1 d Id 2 
cos 9 dr i +sin 9 [rdr + r 2 d0 2 
Id 1 

+2sin0cos<? ( — — 



r 2 d8 r dr dd 



sin 6 d \ 
r dO) 



(8.3) 

(8.4) 



(8.5) 
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and by a similar process, we find that 

3 2 . 2 . 3 2 2 ./l3 13 2 \ 

dy 2 Sm ° dr 2 + C ° S (r dr + r * d ° 2 ) 

-2 S i„»oos e(-L|-iJL) (8.6) 

3 2 / 3 2 13 1 3 2 \ 

“ Sm C ° Se V3r2 “ rdr ~ r^dO 2 J 

-(cos’0- Sin’s) . (8.7) 

Finally, we can determine the expressions for stress components, noting for exam- 
ple that 

a rr = a xx cos 2 0 + a yy sin 2 9 + 2a xy sin 9 cos 9 (8.8) 

, d 2 (b , d 2 <b d 2 <j> 

= cos 2 + sin 2 0— r - 2 sin cos 

ay 2 ox 1 dxdy 

<-> 

after substituting for the partial derivatives from (8.5— 8.7) and simplifying. 

The remaining stress components, ooe, o r e can be obtained by a similar procedure. 
We find 

1 0<j> 1 d 2 <j> d 2 <j> . 

£Jrr “ r 9r + r 2 W ' 006 ~ ( 8-1 °) 



_ 1 d<j> 1 d 2 <j> _ d ( 1 dcj> 

T ° r 2 d6 r dr d9 dr \r dO 



( 8 . 11 ) 



If there is a conservative body force p described by a potential V, the stress 
components are modified to 



13<?J 1 d 2 4> d 2 <(> 

(?rr — 1 2 ^ » a ° e ~ ~a~2 + ^ 

r dr r * 09* dr 1 



- _ 1 
° t8 r 2 d9 r dr 09 



3 / 1 d<j> 
dr 1 r dO 



( 8 . 12 ) 



(8.13) 



Pr = 



Pe = 



We also note that the Laplacian operator 

3 2 3 2 3 2 



V 2 = — ^1-^1 I A iil 

3a; 2 + 3y 2 3r 2 + r 3r + r 2 36 2 



(8.14) 



(8.15) 
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from equations (8.5, 8.6) and hence 

+ ( d 2± + l d A + L?±\ . 

\dr 2 r dr r 2 dO 2 J \dr 2 r dr r 2 dO 2 J 



( 8 . 16 ) 



Notice that in applying the second Laplacian operator in (8.16) it is necessary to 
differentiate by parts. The two differential operators cannot simply be multiplied 
together. In Mathematica and Maple, the easiest technique is to define a new function 



f = = m l_#± 

J ~ 9 dr 2 r dr r 2 dO 2 



and then obtain VV from 



VV = V 2 / 



d ll .!£/ + ±&f 

dr 2 r dr r 2 dO 2 



( 8 . 17 ) 



( 8 . 18 ) 



8.2 Strain components 

A similar technique can be used to obtain the strain-displacement relations in polar 
coordinates. Writing 



u x = u T cos 6 — Uo sin 6 
u y = u T sin 6 + ug cos 6 



and substituting in (8.3), we find 



6tt — 



du x du r 



^ r In, * 

= -7T- COS^ 0 + I 



dx dr 



ug dug 1 du T 



,A 1 dug 

r r dO 



dr 
sin 2 9 



r d6 



sin 6 cos 6 



( 8 . 19 ) 

( 8 . 20 ) 



(8.21) 



Using the same method to obtain expressions for e X y,e yy and substituting the 
results in the strain transformation relations analogous to (8.8) etc., we obtain 



_ du r _ 1 / 1 du r duo uo \ 
dr Cre 2 \r d6 + dr r J 



eee = 



1 due u r 
rd0 + r 



(8.22) 



8.3 Fourier series expansion 

The simplest problems in polar coordinates are those in which there ar e no 0-boundar- 
ies, the most general case being the disk with a central hole illustrated in Figure 8.1 
and defined by b < r < a. 
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Figure 8.1: The disk with a central hole. 



The stresses and displacements must be single-valued and continuous and hence 
they must be periodic functions of 0, since (r, 9 + 2rmr) defines the same point as 
(r, 0), when m is any integer. It is therefore natural to seek a general solution of the 
problem of Figure 8.1 in the form 

00 oo 

<t>=Y, fn( r ) cos (nO) A ffn(f) sin(nfl) . (8.23) 

71=0 71=1 



Substituting this expression into the biharmonic equation, using (8.16), we find 
that the functions /„, g n must satisfy the ordinary differential equation 



d 2 Id n 2 

dr 2 r dr r 2 



(Pfn 1 _ ri 2 f n 

dr 2 r dr r 2 



= 0 , 



which, if n 7^0, 1 has the general solution 

fn(r) = 4ur" +2 + 4 n2 r- n+2 + + >Ur- n , 



(8.24) 



(8.25) 



where A n \, . . . A n 4 are four arbitrary constants. 

Whenn = 0, 1, the solution (8.25) develops repeated roots and equation (8.24) has 
a different form of solution given by 

fo(r) = A 01 r 2 + A 02 r 2 ln(r) + A 03 ln(r) + A 04 (8.26) 

/i( r ) — A\ir 2 A A\ 2 r \n{r) A A\ 3 r A A\ 4 r 1 . (8.27) 

In all the above equations, we note that the stress functions associated with the 
constants A n3 ,A n4 are harmonic and hence biharmonic a fortiori, whereas those as- 
sociated with A n \ , j 4„2 are biharmonic but not harmonic. 
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8.3.1 Satisfaction of boundary conditions 

The boundary conditions on the surfaces r = a,b will generally take the form 

o TT — F x (0) ; r = a (8.28) 

= F 3 (0) ; r = b (8.29) 

v T <) = F 3 (0) ; r = a (8.30) 

= F 4 (B) ; r — b , (8.31) 

each of which has to be satisfied for all values of 9. This can conveniently be done by 
expanding the functions F\, . . . F 4 as Fourier series in 0. i.e. 

oo oo 

Fj(0) = Y a nj cos (n0) + Y b nj sin (n0) ; j = 1, ... ,4 . (8.3 2) 

n=0 n=l 



In combination with the stress function of equation (8.23), (8.28-8.3 T) will then 
give four independent equations for each trigonometric term, cos(n0), sin(n0) in the 
series and hence serve to determine the four constants A nl , . . . , A n4 . The problem of 
Figure 8.1 is therefore susceptible of a general solution. 

8.3.2 Circular hole in a shear field 

The general solution has some anomolous features for the special values n = 0, 1, but 
before discussing these we shall illustrate the method by solving a simple example. 
Figure 8.2 shows a large plate in a state of pure shear o xy = 5, perturbed by a hole 
of radius a. 

s 




Figure 8.2: Circular hole in a shear field. 
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A formal statement of the boundary conditions for this problem is 



tfrr 


= 0 


; r — a 


(8.33) 


&r0 


= 0 


; r = a 


(8.34) 


&xx> ®yy 


-> 0 


; r -> 00 


(8.35) 


Q X y 


-> 5 


; r -A 00 . 


(8.36) 


Notice that we describe a body as 


‘large’ 


when we mean 


to interpret the boundary 



condition at the outer boundary as being applied at r— >oo. Notice also that we have 
stated these ‘infinite’ boundary conditions in (x,y) coordinates, since this is the most 
natural way to describe a state of uniform stress. We shall see that the stress function 
approach makes it possible to use both rectangular and polar coordinates in the same 
problem without any special difficulty. 

This is a typical perturbation problem in which a simple state of stress is perturbed 
by a local geometric feature (in this case a hole). It is reasonable in such cases to 
anticipate that the stresses distant from the hole will be unperturbed and that the 
effect of the hole will only be felt at moderate values of r/a. Perturbation problems 
are most naturally approached by first solving the simpler problem in which the 
perturbation is absent (in this case the plate without a hole) and then seeking a 
corrective solution which will describe the influence of the hole on the stress field. 
With such a formulation, we anticipate that the corrective solution will decay with 
increasing r. 

The unperturbed field is clearly a state of uniform shear, a xy — S, and this in turn 
is conveniently described by the stress function 

<p ~ ~Sxy ~ -Sr 2 sing cos 6 = s ^ 2 — - , (8.37) 

from equation (4.1) 

Notice that although the stress function is originally determined in rectangular 
coordinates, it is easily transformed into polar coordinates using (8.1) and then into 
the Fourier form of equation (8.23). 

The unperturbed solution satisfies the ‘infinite’ boundary conditions (8.35, 8.36), 
but will violate the conditions at the hole surface (8.33, 8.34). However, we can correct 
the stress field by superposing those terms from the series (8.23) which (i) have the 
same Fourier dependence as (8.37) and (ii) lead to stresses which decay as r increases. 
There will always be two such terms for any given Fourier component, permitting the 
two traction boundary conditions to be satisfied. In the present instance, the required 
terms are those derived from the constants 422,^24 in equation (8.25), giving the 
stress function 1 

(j) = - 2~" + 4sin(2<?) + Br~ 2 s\n(26) . (8.38) 



1 It might be thought that the term involving A22 is inappropriate because it does not decay with 
r. However, it leads to stresses which decay with r, which is of course what we require. 
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The corresponding stress components are obtained by substituting in (8.10, 8.11) 
with the result 



o Tr — 


(5 -£-!?)<■(*) 


( 8 . 39 ) 




/ 2 A 6B\ , 




°r9 = 


( s + + -pr) cos ( 20 ) 


( 8 . 40 ) 


°66 = 


(~ 5 + £) sin ( 20 ) • 


( 8 . 41 ) 



The two boundary conditions at the hole surface (8.33, 8.34) then yield the two 
equations 



4A 6B 

~T + — = S 
a 1 or 

2A 6 B 

T + -r = -S 
a 1 a 4 



for the constants A, B, with solution 



A = Sa 2 ■, B = - 



Sa 4 



and the final stress field is 



<r rr = 5 ^1 - 4^ + sin(20) 

o r o = 5^1 + 2^ - 3^ cos(20) 
o eo = 5^-1- 3^-j sin(20) . 



( 8 . 42 ) 

( 8 . 43 ) 

( 8 . 44 ) 

( 8 . 45 ) 

( 8 . 46 ) 

( 8 . 47 ) 



Notice incidentally that the maximum stress is the hoop stress, oee = 45 at the 
point (a, 37t/4). At this point there is a state of uniaxial tension, so the maximum 
shear stress is 25. Since the unperturbed shear stress has a magnitude 5, we say 
that the hole produces a stress concentration of 2. However, for a brittle material, 
we might be inclined to define the stress concentration factor as the ratio of the 
maximum tensile stresses in the perturbed and unperturbed solutions, which in the 
present problem is 4. In general, a stress concentration factor implies a measure of 
the severity of the stress field — usually a failure theory — which serves as a standard 
of comparison and the magnitude of the stress concentration factor will depend upon 
the measure used. 

The determination of the stress concentration due to holes, notches and changes of 
section under various loading conditions is clearly a question of considerable practical 
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importance. An extensive discussion of problems of this kind is given by Savin 2 and 
stress concentration factors for a wide range of geometries are tabulated by Peterson 3 
in a form suitable for use in engineering design. 

8.3.3 Degenerate cases 

We have already remarked in §8.3 above that the solution (8.25) degenerates for 
n = 0, 1 and must be supplemented by some additional terms. In fact, even the 
modified stress function of (8.26, 8.27) is degenerate because the stress function 

(j> = A + Bx + Cy = A + Br cos 0 + Cr sin 0 (8.48) 

defines a trivial null state of stress (see §5.1.1) and hence the constants A 0 4 ,Ai 3 in 
equations (8.26, 8.27) correspond to null stress fields. 

This question of degeneracy arises elsewhere in Elasticity and indeed in mathe- 
matics generally, so we shall take this opportunity to develop a general technique for 
resolving it. As we saw in §8.3, the degeneracy often arises from the occurrence of 
repeated roots to an equation. Thus, in equation (8.25), the terms A„ 2 r~" +2 , A n3 r n 
degenerate to the same form when n = 1 . 

Suppose for the moment that we relax the restriction that n be an integer. Clearly 
the degeneracy in equation (8.25) only arises exactly at the values n = 0, 1. There is 
no degeneracy for n= 1+e for any non-zero e, however small. We shall show therefore 
that we can recover the extra solution at the degenerate point by allowing the solution 
to tend smoothly to the limit n -4 I , rather than setting n = 1 ab initio. 

For ?r=l+e, the two offending terms in (8.25) can be written 

f(r) = Ar l ~ l + Dr 1+e , (8.49) 

where A. B, are two arbitrary constants. 

Clearly the two terms tend to the same form as e -4 0. However, suppose we 
construct a new function from the sum and difference of these functions in the form 

f{r) = C(r l+C + r'~ c ) + D{r l+e - r l ~‘) . (8.50) 

In this form, the first function tends to 2 Cr as e -4 0, whilst the second tends 
to zero. However, we can prevent the second term degenerating to zero, since the 
constant D is arbitrary. We can therefore choose D = Et~ l in which case the second 
term will tend to the limit 



.Hi _ ~i-< 

lim E = 2£rln(r) , 

<— >o e 



(8.51) 



2 G.N. Savin, St?'ess Concentration around Holes, Pergamon Press, Oxford, (1961). 

3 R.E. Peterson, Stress Concentration Design Factors, John Wiley, New York, (1974). 
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where we have used L’HopitaFs rule to evaluate the limit. Notice that this result 
agrees with the special term included in equation (8.27) above to resolve the degen- 
eracy in <p for n = 1 . 

This procedure can be generalized as follows: Whenever a degeneracy occurs at a 
denumerable set of values of a parameter (such as n in equation (8.25)), we can always 
make up the defecit using additional terms obtained by differentiating the original 
form with respect to the parameter before allowing it to take the special value. 

The pairs r n+2 ,r~ n+2 and " n , r~ n both degenerate to the same form for n = 0, so 
the new functions are of the form 

rl(r n+2 ) 

la-V = rMn(r > (8 ' 52) 

Haisr = ln(r) - (853) 

again agreeing with the special terms introduced in (8.26). 

Sometimes the degeneracy is of a higher order, corresponding to more than two 
identical roots, in which case L’Hopital’s rule has to be applied more than once — i.e. 
we have to differentiate more than once with respect to the parameter. It is always a 
straightforward matter to check the resulting solutions to make sure they are of the 
required form. 

We now turn our attention to the degeneracy implied by the triviality of the stress 
function of equation (8.48). Here, the stress function itself is not degenerate — i.e. it 
is a legitimate function of the required Fourier form and it is linearly independent of 
the other functions of the same form. The trouble is that it gives a null stress field. 
We must therefore look for stress functions that are not of the standard Fourier form, 
but which give Fourier type stress components. 

Once again, the method is to approach the solution as a limit using L’Hopital’s 
rule, but this time we have to operate on the complete stress function — not just on 
the part that varies with r. Suppose we consider the axisymmetric degenerate term, 
(j)=A, which can be regarded as the limit of a stress function of the form 

0 = Ar e cos(€0) + Br e sin(f0) , (8.54) 

as e-»0. 

Differentiating this function with respect to e and then letting t tend to zero, we 
obtain the new function 

<$> = Aln(r) + B6 . (8.55) 

The first term is the same one that we found by operating on the function / n (r) 
and is of the correct Fourier form, but the second term is not appropriate to a Fourier 
series. However, when we substitute the second term into equations (8.10, 8.11), we 
obtain the stress components 

B 

O r r — Ct 00 — 0 j &r8 — > 

r 2 



(8.56) 




106 



CHAPTER 8. PROBLEMS IN POLAR COORDINATES 



which are of the required Fourier form (for n — 0), since the stress components do not 
vary with 9. 

In the same way, we can develop special terms to make up the deficit due to the 
two null terms with n=l (equation (8.48), obtaining the new stress functions 

$ = Br6 sin 6 + Cr9 cos 9 , (8.57) 



which generate the Fourier-type stress components 

IB cos 9 2Csin0 

@r0 ~ @00 ~ 0 , @rr ~ 

r r 



(8.58) 



8.4 The Michell solution 

The preceding results now permit us to write down a general solution of the elasticity 
problem in polar coordinates, such that the stress components form a Fourier series 
in 0. We have 

<p = -(- A 02 r 2 ln(r) -I- A 0 3 ln(r) + A Q4 9 

+(A u r 3 + A\ 2 r ln(r) + Ai 4 r _1 ) cos 0 + yli 3 r0sin 0 
+(Bur 3 + B\ 2 r ln(r) + 73i 4 r _l ) sin 9 + B i3 r9 cos 9 

OO 

+ ZK j4 «i r ' l+2 + 7l n2 r~ n+2 + A n3 r n + A n4 r~ n ) cos (nO) 

n-2 

+ f^(B nl r n+2 + B n2 r~ n+2 + B n3 r n + B n ,r~ n ) sin (n9) (8.59) 

n=2 

This solution is due to Michell 4 . The corresponding stress components are easily 
obtained by substituting into equations (8.10, 8.11). For convenience, we give them 
in tabular form in Table 8.1, since we shall often wish to select a few components from 
the general solution in the solution of specific problems. The terms in the Michell 
solution are also given in the Maple and Mathematica files ‘Michell’, from which 
appropriate terms can be cut and pasted as required. 

Notice that there are four independent stress functions for each term in the Fourier 
series, as required by the argument of §8.3.1. If the disk is solid, there is no inner 
boundary and we must exclude those components in Table 8.1 that give stresses which 
go to infinity as r — > 0. 

In addition to the special stress functions necessitated by the degeneracy discussed 
in §8.3.3, the cases n = 0, 1 exhibit other anomolies. In particular, some of the stress 
functions for n = 0, 1 correspond to multiple-valued displacement fields and cannot 
be used for the complete annulus of Figure 8.1. Also, equilibrium requirements place 

4 J.H.Michell, On the direct determination of stress in an elastic solid, with application to the 
theory of plates, Proc, London Math. Soc., Vol. 31 (1899), 100-124. 
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Table 8.1: The Michell solution — stress components 



<f> 


<7 rr 


OtO 


<7 00 


r 2 


2 


0 


2 


r 2 ln(r) 


2 ln(r) + 1 


0 


2 ln(r) + 3 


ln(r) 


1/r 2 


0 


— 1/r 2 


# 


0 


1/r 2 


0 


r 3 cos 9 


2 r cos 9 


2r sin# 


6r cos # 


r9 sin 9 


2 cos 9/r 


0 


0 


r ln(r) cos 9 


cos 9/r 


sin 9/r 


cos 9/r 


cos 9/r 


-2 cos 9/r 3 


-2 sin 9/r 3 


2 cos 9/r 3 


r 3 sin 9 


2r sin 9 


—2 r cos# 


6r sin # 


r9 cos 9 


2 sin 9/r 


0 


0 


r ln(r) sin# 


sin 9/r 


— cos 9/r 


sin 9/r 


sin#/r 


-2 sin 9/r 3 


2 cos 9/r 3 


2 sin 9/r 3 


r ,H2 cos n# 


— (n+l)(n— 2)r n cos n9 


n(n+l)r"sin n# 


(n + 1 ) (n+ 2)r"cos n# 


r -«+ 2 cos n 9 


— (n+2)(n— l)r -n cosn# 


— n(n— l)r _n sin n# 


(n — l)(n — 2)r _n cos n# 


r n cos n9 


-n(n— l)r"- 2 cos n9 


n(n— ljr^sin n# 


n(n— ljr^cosn# 


r ~ n cos n9 


— n(n+ l)r _n_2 cos n9 


— n(n4-l)r _n_2 sin n# 


n(n+l)r _n_2 cosn# 


r^sin n9 


— (n+l)(n — 2)r"sin n9 


— n(n+l)r n cosn# 


(n+l)(n+2)r"sin n# 


r _Itf2 sin n9 


— (n + 2)(n— l)r - "sin n9 


n(n-l)r“”cosn# 


(n— l)(n— 2)r _n sin n# 


r"sin n9 


— n(n— l)r"^sin n9 


— n(n— l)r"" 2 cos n# 


n(n— l)r" _2 sin n9 


r _n sin ti9 


— n(n+l)r _n_2 sin n9 


n(n+l)r _,l_2 cos n# 


n(n + l)r _n_2 sin n# 
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restrictions on the permissible boundary conditions for the terms n = 0, 1 in the 
Fourier series (8.32). 

Detailed discussion of these difficulties will be postponed to §9.3.1, after we have 
introduced methods of determining the displacements associated with a given stress 
field. In the present chapter, we shall restrict attention to cases where these difficulties 
do not arise. 

8.4.1 Hole in a tensile field 

To illustrate the use of Table 8.1, we consider the case where the the body of Figure 
8.2 is subjected to uniform tension at infinity instead of shear, so that the boundary 
conditions become 



o rr 


— 


o ; 


r = a 


(8.60) 


O r Q 


= 


0 ; 


r = a 


(8.61) 


a w 




o ; 


r oo 


(8.62) 
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(8.63) 



The unperturbed problem in this case can clearly be described by the stress func- 
tion 



. Sy 2 St 2 sin 2 9 Sr 2 Sr ! cos(20) 

— = — 2 4 4 



(8.64) 



This function contains both an axisymmetric term and a cos(20) term, so to 
complete the solution, we supplement it with those terms from Table 8.1 which have 
the same form and for which the stresses decay as r-4oo. 

The resulting stress function is 




Sr 2 cos(20) 
4 



+ A ln(r) + B6 + C cos 26 + Dr 2 cos(2 6) 



(8.65) 



and the corresponding stress components are 

S S cos(20) A 4Ccos(20) 6Dcos(20) 

= 2 + — + ^ — ? ;; — (8 - 66) 

Ssin(20) B 2Csin(20) 6Dsin(20) 

°re = j + ^ r* (8 ’ 67) 

S S cos(20) A 6D cos(26>) 

^ 2 ^ + r* ' (868) 



The boundary conditions (8.60, 8.61) will be satisfied if and only if the coefficients 
of both Fourier terms are zero onr = o and hence we obtain the equations 




0 



(8.69) 




